{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 线性回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 初始化模型参数\n",
    "def initialize_params(dims):\n",
    "    '''\n",
    "    输入：\n",
    "    dims：训练数据变量维度\n",
    "    输出：\n",
    "    w：初始化权重参数值\n",
    "    b：初始化偏差参数值\n",
    "    '''\n",
    "    # 初始化权重参数为零矩阵\n",
    "    w = np.zeros((dims, 1))\n",
    "    # 初始化偏差参数为零\n",
    "    b = 0\n",
    "    return w, b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "### 定义模型主体部分\n",
    "### 包括线性回归公式、均方损失和参数偏导三部分\n",
    "def linear_loss(X, y, w, b):\n",
    "    '''\n",
    "    输入:\n",
    "    X：输入变量矩阵\n",
    "    y：输出标签向量\n",
    "    w：变量参数权重矩阵\n",
    "    b：偏差项\n",
    "    输出：\n",
    "    y_hat：线性模型预测输出\n",
    "    loss：均方损失值\n",
    "    dw：权重参数一阶偏导\n",
    "    db：偏差项一阶偏导\n",
    "    '''\n",
    "    # 训练样本数量\n",
    "    num_train = X.shape[0]\n",
    "    # 训练特征数量\n",
    "    num_feature = X.shape[1]\n",
    "    # 线性回归预测输出\n",
    "    y_hat = np.dot(X, w) + b\n",
    "    # 计算预测输出与实际标签之间的均方损失\n",
    "    loss = np.sum((y_hat-y)**2)/num_train\n",
    "    # 基于均方损失对权重参数的一阶偏导数\n",
    "    dw = np.dot(X.T, (y_hat-y)) /num_train\n",
    "    # 基于均方损失对偏差项的一阶偏导数\n",
    "    db = np.sum((y_hat-y)) /num_train\n",
    "    return y_hat, loss, dw, db"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 定义线性回归模型训练过程\n",
    "def linear_train(X, y, learning_rate=0.01, epochs=10000):\n",
    "    '''\n",
    "    输入：\n",
    "    X：输入变量矩阵\n",
    "    y：输出标签向量\n",
    "    learning_rate：学习率\n",
    "    epochs：训练迭代次数\n",
    "    输出：\n",
    "    loss_his：每次迭代的均方损失\n",
    "    params：优化后的参数字典\n",
    "    grads：优化后的参数梯度字典\n",
    "    '''\n",
    "    # 记录训练损失的空列表\n",
    "    loss_his = []\n",
    "    # 初始化模型参数\n",
    "    w, b = initialize_params(X.shape[1])\n",
    "    # 迭代训练\n",
    "    for i in range(1, epochs):\n",
    "        # 计算当前迭代的预测值、损失和梯度\n",
    "        y_hat, loss, dw, db = linear_loss(X, y, w, b)\n",
    "        # 基于梯度下降的参数更新\n",
    "        w += -learning_rate * dw\n",
    "        b += -learning_rate * db\n",
    "        # 记录当前迭代的损失\n",
    "        loss_his.append(loss)\n",
    "        # 每1000次迭代打印当前损失信息\n",
    "        if i % 10000 == 0:\n",
    "            print('epoch %d loss %f' % (i, loss))\n",
    "        # 将当前迭代步优化后的参数保存到字典\n",
    "        params = {\n",
    "            'w': w,\n",
    "            'b': b\n",
    "        }\n",
    "        # 将当前迭代步的梯度保存到字典\n",
    "        grads = {\n",
    "            'dw': dw,\n",
    "            'db': db\n",
    "        }     \n",
    "    return loss_his, params, grads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(442, 10)\n",
      "(442,)\n",
      "[[ 0.03807591  0.05068012  0.06169621  0.02187235 -0.0442235  -0.03482076\n",
      "  -0.04340085 -0.00259226  0.01990842 -0.01764613]\n",
      " [-0.00188202 -0.04464164 -0.05147406 -0.02632783 -0.00844872 -0.01916334\n",
      "   0.07441156 -0.03949338 -0.06832974 -0.09220405]\n",
      " [ 0.08529891  0.05068012  0.04445121 -0.00567061 -0.04559945 -0.03419447\n",
      "  -0.03235593 -0.00259226  0.00286377 -0.02593034]\n",
      " [-0.08906294 -0.04464164 -0.01159501 -0.03665645  0.01219057  0.02499059\n",
      "  -0.03603757  0.03430886  0.02269202 -0.00936191]\n",
      " [ 0.00538306 -0.04464164 -0.03638469  0.02187235  0.00393485  0.01559614\n",
      "   0.00814208 -0.00259226 -0.03199144 -0.04664087]]\n",
      "[151.  75. 141. 206. 135.]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.datasets import load_diabetes\n",
    "diabetes = load_diabetes()\n",
    "data = diabetes.data\n",
    "target = diabetes.target \n",
    "print(data.shape)\n",
    "print(target.shape)\n",
    "print(data[:5])\n",
    "print(target[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X_train's shape:  (353, 10)\n",
      "X_test's shape:  (89, 10)\n",
      "y_train's shape:  (353, 1)\n",
      "y_test's shape:  (89, 1)\n"
     ]
    }
   ],
   "source": [
    "# 导入sklearn diabetes数据接口\n",
    "from sklearn.datasets import load_diabetes\n",
    "# 导入sklearn打乱数据函数\n",
    "from sklearn.utils import shuffle\n",
    "# 获取diabetes数据集\n",
    "diabetes = load_diabetes()\n",
    "# 获取输入和标签\n",
    "data, target = diabetes.data, diabetes.target \n",
    "# 打乱数据集\n",
    "X, y = shuffle(data, target, random_state=13)\n",
    "# 按照8/2划分训练集和测试集\n",
    "offset = int(X.shape[0] * 0.8)\n",
    "# 训练集\n",
    "X_train, y_train = X[:offset], y[:offset]\n",
    "# 测试集\n",
    "X_test, y_test = X[offset:], y[offset:]\n",
    "# 将训练集改为列向量的形式\n",
    "y_train = y_train.reshape((-1,1))\n",
    "# 将验证集改为列向量的形式\n",
    "y_test = y_test.reshape((-1,1))\n",
    "# 打印训练集和测试集维度\n",
    "print(\"X_train's shape: \", X_train.shape)\n",
    "print(\"X_test's shape: \", X_test.shape)\n",
    "print(\"y_train's shape: \", y_train.shape)\n",
    "print(\"y_test's shape: \", y_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch 10000 loss 3679.868273\n",
      "epoch 20000 loss 3219.164522\n",
      "epoch 30000 loss 3040.820279\n",
      "epoch 40000 loss 2944.936608\n",
      "epoch 50000 loss 2885.991571\n",
      "epoch 60000 loss 2848.051813\n",
      "epoch 70000 loss 2823.157085\n",
      "epoch 80000 loss 2806.627821\n",
      "epoch 90000 loss 2795.546917\n",
      "epoch 100000 loss 2788.051561\n",
      "epoch 110000 loss 2782.935842\n",
      "epoch 120000 loss 2779.411265\n",
      "epoch 130000 loss 2776.957989\n",
      "epoch 140000 loss 2775.230803\n",
      "epoch 150000 loss 2773.998942\n",
      "epoch 160000 loss 2773.107192\n",
      "epoch 170000 loss 2772.450534\n",
      "epoch 180000 loss 2771.957489\n",
      "epoch 190000 loss 2771.579121\n",
      "{'w': array([[  10.56390075],\n",
      "       [-236.41625133],\n",
      "       [ 481.50915635],\n",
      "       [ 294.47043558],\n",
      "       [ -60.99362023],\n",
      "       [-110.54181897],\n",
      "       [-206.44046579],\n",
      "       [ 163.23511378],\n",
      "       [ 409.28971463],\n",
      "       [  65.73254667]]), 'b': 150.8144748910088}\n"
     ]
    }
   ],
   "source": [
    "# 线性回归模型训练\n",
    "loss_his, params, grads = linear_train(X_train, y_train, 0.01, 200000)\n",
    "# 打印训练后得到模型参数\n",
    "print(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 82.0537503 ],\n",
       "       [167.22420149],\n",
       "       [112.38335719],\n",
       "       [138.15504748],\n",
       "       [174.71840809]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### 定义线性回归预测函数\n",
    "def predict(X, params):\n",
    "    '''\n",
    "    输入：\n",
    "    X：测试数据集\n",
    "    params：模型训练参数\n",
    "    输出：\n",
    "    y_pred：模型预测结果\n",
    "    '''\n",
    "    # 获取模型参数\n",
    "    w = params['w']\n",
    "    b = params['b']\n",
    "    # 预测\n",
    "    y_pred = np.dot(X, w) + b\n",
    "    return y_pred\n",
    "# 基于测试集的预测\n",
    "y_pred = predict(X_test, params)\n",
    "# 打印前五个预测值\n",
    "y_pred[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 37.]\n",
      " [122.]\n",
      " [ 88.]\n",
      " [214.]\n",
      " [262.]]\n"
     ]
    }
   ],
   "source": [
    "print(y_test[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 定义R2系数函数\n",
    "def r2_score(y_test, y_pred):\n",
    "    '''\n",
    "    输入：\n",
    "    y_test：测试集标签值\n",
    "    y_pred：测试集预测值\n",
    "    输出：\n",
    "    r2：R2系数\n",
    "    '''\n",
    "    # 测试标签均值\n",
    "    y_avg = np.mean(y_test)\n",
    "    # 总离差平方和\n",
    "    ss_tot = np.sum((y_test - y_avg)**2)\n",
    "    # 残差平方和\n",
    "    ss_res = np.sum((y_test - y_pred)**2)\n",
    "    # R2计算\n",
    "    r2 = 1 - (ss_res/ss_tot)\n",
    "    return r2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5334188457463576\n"
     ]
    }
   ],
   "source": [
    "print(r2_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAELCAYAAAAspXpuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsvXmYHFd57/95Z99nNJrRNrIs2Zbl3ZItjDfAGO8sFgQCBAIhBBNitoT4YpPf74bkXgcHJ+wJiVkNgbAFCwMGAzYYMNhYtmRbi2Vrl0bbzEiz7zPn/nHqdFdXV1VXb9Pd0+fzPPP0dE11d0111XnP+32XI0opLBaLxWLxo6LQB2CxWCyW4sUaCYvFYrEEYo2ExWKxWAKxRsJisVgsgVgjYbFYLJZArJGwWCwWSyDWSFgsFoslEGskLBaLxRKINRIWi8ViCaSq0AeQLR0dHWrlypWFPgyLxWIpGZ588slepVRnlH1L3kisXLmSTZs2FfowLBaLpWQQkf1R97Vyk8VisVgCsUbCYrFYLIHkzUiISJ2I/EFEnhaRbSLyD872r4rIXhHZ4vysdbaLiHxGRHaJyDMiclG+js1isVgs0chnTGICuFopNSwi1cBvReQnzt9uU0p9z7P/jcBq5+fFwOedR4vFYrEUiLx5Ekoz7Dytdn7CFq+4Gfia87rHgDYRWZqv47NYLBZLavIakxCRShHZAhwHfq6Uetz5052OpPRJEal1tnUBB10vP+RssxQxGzd3c8VdD7Pq9h9zxV0Ps3Fzd6EPyWKx5JC8Ggml1IxSai2wHLhERM4D7gDOAl4EtAMfdnYXv7fwe18RuUVENonIpp6enjwcuSUKGzd3c8f3n6W7fwwFdPePccf3n7WGwmKZR8xJdpNSqh/4FXCDUuqIIylNAF8BLnF2OwSc4nrZcuBwwPvdo5Rar5Ra39kZqR7EkgfufnAnY1MzCdvGpma4+8GdBToii8WSa/KZ3dQpIm3O7/XANcBzJs4gIgJsALY6L7kfeJuT5XQpMKCUOpKv47Nkz+H+sbS2WyyW0iOf2U1LgXtFpBJtjL6jlPqRiDwsIp1oeWkL8JfO/g8ANwG7gFHgHXk8NksOWNZWT7ePQVjWVl+Ao7FYLPkgb0ZCKfUMsM5n+9UB+yvg1nwdjyX33Hb9Gu74/rMJklN9dSW3Xb+mgEdlsVhyScn3brIUjg3rdPLZ3Q/u5HD/GMva6rnt+jWx7RaLpfSxRsKSFRvWdRXMKGzc3G0NlMWSZ6yRsJQkJv3WSF0m/RawhsJiySHWSFhKkrD0W2sk/LGelyUTrJGwlCQ2/TY9rOdlyRTbKtxSkgSl2dr0W39s4aMlU6yRsJQkt12/hvrqyoRtNv02GOt5WTLFGglLSbJhXRcfe935dLXVI0BXWz0fe935VjoJwHpelkyxMQlLyVLI9NtSwxY+WjLFGgmLpQywhY+WTLFGwmIpE6znZckEG5OwWCwWSyDWkyhSbOGTxWIpBqyRKEJs4ZPFYikWrNxUhNjCJ4vFUixYI1GE2MIni8VSLFgjUYTYwieLxVIsWCNRhNiWExaLpViwgesixBY+FQ6bVWaxJGKNRJFiC5/mHptVZrEkY+Umi8XBZpVZLMlYI2GxONisMoslmbwaCRGpE5E/iMjTIrJNRP7B2b5KRB4XkRdE5NsiUuNsr3We73L+vjKfx2exuLFZZRZLMvn2JCaAq5VSFwJrgRtE5FLgn4FPKqVWAyeBdzr7vxM4qZQ6A/iks5/FMifYrDKLJZm8GgmlGXaeVjs/Crga+J6z/V5gg/P7zc5znL+/QkQkn8dosRjsQkYWSzJ5z24SkUrgSeAM4N+A3UC/Umra2eUQYO7CLuAggFJqWkQGgIVAb76P0xJMOaWFZppVVk7nyFJe5N1IKKVmgLUi0gbcB5ztt5vz6Oc1KO8GEbkFuAVgxYoVOTpSix82LTQ19hxZ5jNzlt2klOoHfgVcCrSJiDFQy4HDzu+HgFMAnL+3Aid83usepdR6pdT6zs7OfB96WWPTQlNjz5FlPpPv7KZOx4NAROqBa4AdwC+B1zu7vR34gfP7/c5znL8/rJRK8iQsc0e+0kI3bu7mirseZtXtP+aKux5m4+burN6vkNjUWct8Jt9y01LgXicuUQF8Ryn1IxHZDnxLRP4vsBn4krP/l4Cvi8gutAfxpjwfnyUFy9rq6fYZ7LJJC51v8kw+zpHFUizkO7vpGaXUOqXUBUqp85RS/+hs36OUukQpdYZS6g1KqQln+7jz/Azn73vyeXyW1OQjLXS+yTM2ddYyn7G9myyh5KPZ4HyTZ2xDRst8xhoJS0py3WxwPsoztiGjZb5iezdZ5hwrz1gspYP1JCxzjpVnLG5sIWJxY42EpSBYecYC8y/TbT5i5SaLxVIw5lum23zEGgmLxZLMwF6471UwNZLXj5lvmW7zESs3WSyWZA79Bvb8GE7ugkUX5u1j0sl0s7GLwmA9iXnOfGp/YZlDJof048x4Xj8maqabiV1094+hiMcu7PWcf6wnMY+xQUFLxsSMxERePyZqpltY7KKUruVS9IaskZjHzJcbyzK3bNzczcAjz/L2SvjgNx7jqutOz+v1EiXTbT7ELkp10mblpnnMfLixLHOLGchmJ7QnMTQyXBSyznxYf7xUM7mskZjHzIcbyzK3mIGsuWIUgBqmimIgmw9V+qU6abNGYh4zH24sy9xiBqxG0Y81MpWwPYx8JknMh/XHS3XSZmMS8xjb/sKSLiYl1RiJWsdIpBrI5kJvL/Uq/duuX5NwjqA0Jm3WSKRBKWYmlPqNZZlbzEDWJHG5KcpAZpMkUlOqkzZrJCJSqpkJFks6mGu57SGd+rqoQfGxa4NlHTNx8iuIg+LX2+eaUpy02ZhEREo1M8FiSZcN67o4rWUWgPe/bEWogTAFbkEUu95uSY01EhEp1cwEiyUjplIX0/lNnNyUgt5uSY01EhEp1cwEiyVtlIpXXE8Ht+UImyCVYvaRxR9rJCJi00ktZcPMBMxOx38PIGiC1NVWz6O3X20NxDzBGomIzIc87bxx/+vh6f8M3cU2GiwhjBcBoUbCTpzKg7xlN4nIKcDXgCXALHCPUurTIvJR4F1Aj7PrR5RSDzivuQN4JzADvF8p9WC+ji8TSjEzYU7Y91OoboQL3+37Z5sZVmK4jUSI3FSqKZ2W9MhnCuw08CGl1FMi0gw8KSI/d/72SaXUv7h3FpFzgDcB5wLLgF+IyJlKqeDImKXwzEzphWmmRwN3sTn0JUZETwLsxKkcyJvcpJQ6opR6yvl9CNgBhF1NNwPfUkpNKKX2AruAS/J1fJYcMTGgH0NWMCu5zLBd98NoT+r95iuTg/Hf89wq3FL8zElMQkRWAuuAx51N7xWRZ0TkyyKywNnWBRx0vewQAUZFRG4RkU0isqmnp4xv5mJgol8/TgV7EiWVGTY5BD/YAM9+odBHUjgSPIn8LjpkKX7ybiREpAn4H+CDSqlB4PPA6cBa4Ajwr2ZXn5crv/dUSt2jlFqvlFrf2dmZh6O2RCZmJII9iZIKcI72AApGjxf0MAoa6DdGorYVpq0nUe7ktS2HiFSjDcQ3lFLfB1BKHXP9/QvAj5ynh4BTXC9fDhzO5/EVmo2bu/m3Bzfz8on7+GntG/mb688tPX3XGImQmERJBTjH+/TjWG/BDqHggX5jJOo7rCdhyWt2kwBfAnYopT7h2r5UKXXEefpaYKvz+/3AN0XkE+jA9WrgD/k6vkJjBoKXy6N8pPXL7BpYzh3f160QinLwDCKCJwElFOAcK7yRKHigf8ptJKwnUe7kU266AvhT4GoR2eL83AR8XESeFZFngJcDfw2glNoGfAfYDvwUuHU+ZzbFF3fRg+uV1VtKsxfUeOqYRElRBJ5EwQP9kwUwEiPHoH/P3HyWJS3y5kkopX6Lf5zhgZDX3Ancma9jKibMDW9aMl9RswVGijjjJ4iInkTe+Pm7oaYFXnZ3bt4voieRz7bxZk0Hv+1zwuQQVNVDVWNonUROeeRvofcZeNvT4fsN7ocH/xxe/T2oWxC+7zyikMsU2IrrAmFu+GbHSKypOkBnxQkqREqrKtkdk1Czc//5Bx+B/T/L3ftFMBLu7qeKeMwgV99XwQP9k0NQ0wxVdXPnSYweg/7dum9UGN2PwoGHoW/73BxXEZDv6y0V1kgUCDMQGE8CtOQ0o1RBLoSMMUYC5m7W6WasV88uc4WRm6ZGYMrfq8t32/iCt4AxRqKydu6MxOSwPufuGg0/TNbZdIl53FlQ6GUK7KJDBcLc8JW/mOTYTDtVMs2VNVu4b+Lq2D4lUZXsNhJTI1DdMHefPTsD4ycApYv6aluzf0/jSYA2GNXLk3aZi5hBQQP9k0NQ7RiJuTL8U8P6cbibjduHg6WVMacuqoyMRKFjVNaTKCAb1nXx6rMaWdy5mN9NXsgV1VvwloYUfYwiwZOY4+C1MRCQO2/CLTMFSE4lVRyYCVMFkJscI/Ho5i3h0oqphJ8viRIRKPT1Zo1EoZkYhNoWtlW9iCWVJzi98lDCn4t+4Bn3eBJzyZir2j5XRmK8Tw+QEGgkCh4zyDdeuSlVnCBXnwn86skt4dJKGcpNhb7erJEoNJODUNPCRS95A6DjEoacXAhKwaZP5K8X0UQ/VFTr3+d6ducexHPmSfTBgjXJ7++i4DGDfDM5GPck1Gx8bYl84ngSdRNHfP8c86jLUG4q9PVmYxKFZnIQmrq47vJLGd5yKtfMPsvXxl+duzS3wX3wyId0SuPa9+TkkBOY6IemZXqQnnNPIg9GYrwPlr8Ujm0KzXAqmeLATHB7EqC9icrq/H3e7HQs9nF6XT/4zDViHnXMkygfuQkKe71ZI1FoJrQnAdC0+gZe8tw32ftP10NFjr6aCSdbxHRrzTUT/dB5oR6k5/rGNd5RTbM2htkyPaEN3YLV+nkBC+oKitdITI9DTVMeP2849uuli8apH6pMkJwSPOoy9CQKjZWbCs2kjkkAcOo1+gY9/Bjse1Cv+Pa1ddnJOCalcDIPRsKsJdG0TD8vlCex6KLceBIm/bVhsS7Ucmc6lQtq1slSc+QmyH/weipuJJZU9gVLKzOT8cmONRJzhvUkComadWZtTurmKS8HBP7nWj17q6zRN8aJnbB4XWafYVosTKTIP88Ec8M2OW5wIWISNc165r/7/hy8n2MU6hfqlhTl6EmYWb1Xbgqjb7u+BjJNQXZ/5tChYGnFHVcro+ymQmM9iUIyOQyouCdRvxDOfRt0vQRe+S3441/p7dlIKcZI5MOTMOmvjYXyJHr0YN5yqtaqA4rfImM8ibqFUFeuRsK0CW9xGYkUtRLfuhKe+Hjmn2kaCi5Yo7/ToPbk7my2eeBJlMq679aTKCRGCnJiEgDc8NX472Mn9OPA3uw/Ix+ehDE8zcaTKIDcVN+pjQTA0AFozyIbzBgF40kMHQzffz5ijER1M1Q6clPYmhIzUzB+EgYPZPGZjifRfpZOGBg5Aq0rk/dzr/FR4kai4O3g08B6EoXEz0i4qVug/1asnsS4x5PwC1xv/jfo25H7zwYtPxhPArKPS4y5PIlylZvMrL6mGaoiyE3m+spmkaYpl5EAGA6YURtPoqKq5LObCt1qIx2skSgkZnZfG2AkRPSMamBf5p+R15iEYyTqF2ppwutJzE7Dw++FZ7+U+88Gx5PIoZEY98QkxnvnppCsmJh0GYkocpMxKmNZ1OGYzzRe4NAh//1MTKJpecl7EoVutZEOVm4KIe/teVN5EgAtK3MjN+UzJlHbBtWNycFEc/OP5ylLyMhNTV0glbnxJKoadFZPfYdOHpge1f+bH/se1Ab8wndn97nFhNtImL5NxeJJjB7XXkTj0pI3Etm0g5/rtuHWkwhgTtrzRjUSg/syn9HOhSdR26YHV68nYbKf8iHbTI3qAby+Qw8cTV258STqF+rf6zv0Y9ixP/2f8PD7YLSEZKnZmfCW7pM+clNYTMJtJDK+Rh0j0XyKLvoMk5vqO/wnJCVGpq02CtE23BqJAOZEM0wlNwG0rtI34viJzD7DHZPIRjrZ8U3Y51m3YaIfpEIXWlU3JuvExgjmw0iY92zo1I+tK3PjSRjjYIxF2LFPj8LsFDz3jew+dy75+S3wg9cG/z02cXEFrsPkJrP/7FTmBZuxYHmTNvaBctNxaFikDUmJexKZttooRCzDyk0BzIlmaCSgmpD88paV+nFwX3zgSgejGZvWB9UZNgx87P/oG3jldfFt4/06N14qdIvwIE8iH3JTLBPJGdRbTtULEGX1nn06aO1+31Aj4VwLW78CF30gu8+eKw7/HsJWBXZ7EmaGH0VuAj2I17Wlf0xTw3rgr6jU11ig3NSj5cXqhoIbiVxIPpm02ihELMN6EgHMSXveCdesLQiTCphp8NotM2UTl5gc1immCe/dr6UmKIAn4QQx3UZi+JBOycyU8d705KapUUCg52k4tjnzz50r1CwM7Amf8U8OaaNf1ZDYliNsf0OmcYmp4fg90Lw8hdzUqQ1KAeWmQq4UV4i24dZIBDAn7XknB/XgWlEZvI/bk8joM1w3cTZxialhnQvv1rPdRsI3JuF83vjJjDuJBhYcxTwJR25qPlUfW9AAE4W0PYlROOUqXRm/7SuZf+5cMXxYewVhq79NDmnZRyRaWw739ZVphtPksP5MiHsSfnGTIpGbCpm+Woi24dZIBDAn7XknB8OD1qDd99q2zDOcpobin5GpJ6GUNhIzE4mtEbyehNdIuD8vg5hK6IzNT26CzOMSszPamBlPorZNz6hTyU3Ny+GM18KOb4QHeIuB/t36cXos2OMyzf0gWluOXHgSk0PxBoJNXTq+4T3v045xa+jUE5ICGolCpq8Wom14ypiEiHxdKfWnqbb5vO4U4GvAEmAWuEcp9WkRaQe+DawE9gF/rJQ6KSICfBq4Cd0s+M+UUk+l/y/ljny053VrmV9cuItLmhoIEZs0JsMpEyYG9UDWtz3zwOLMZNwTGNwPjYud9+6HNqdjql/GidtzGevTs8A0CJuxbbiyR6e9Gg08WyMx0Q+ouJGoqIS69tRyU1U9nP0W2Plt2PNDOPP1mX3+XGCMBOhroaEjeR/3pCKS3DSos8tmp7OTm6pdchPAUHfi9RKTFzu10Zoe1ZMXkcw+MwuySV/NBXPdNjyKJ3Gu+4mIVAIXR3jdNPAhpdTZwKXArSJyDnA78JBSajXwkPMc4EZgtfNzC/D5SP9BCeGdGVdOD7F3sDK1lplNQd3UkC4+gtSLzAfhni264xLjbk/CR25yexIZxCVCZ2xjTvxAnEu4ZYV+zNRIuKutDamqrqdH9ax2xTX6HG/NTnKK2ssn454//bvivwd5lQmeRI1+TOVJ1LTq6yCrmITLkwAdX3JjPFgjN6U6rjxS6JXi5ppAIyEid4jIEHCBiAw6P0PAceAHqd5YKXXEeAJKqSFgB9AF3Azc6+x2L7DB+f1m4GtK8xjQJiJLM/3HihHvzLhJRhmcrU+tZWZaKzEzpWeB5sbL1JNwtXJOGIQn+uMzeb/AdYInkb6RCA3SmZx5Q1WdbvGdppEwA+7rPqG7yP7usOsc13eEtwufHtPGsaJSN2bc91Ot+2dA1GBoVkFTryfhh9tIiMSXMA1ickincDcsyk5uqvYaCc//4/Ykqhr07xElp1w30iv0SnFzTaCRUEp9TCnVDNytlGpxfpqVUguVUnek8yEishJYBzwOLFZKHXE+4whgfMouwN1R7ZCzbd7gnRk3ywhDqjG1ltm6Us/S0x1ojQfQnKUnkWAkHE9idlpv9wau3YZscjC+tGkGRiJ0xmaqrd20nJqWkXAPuG0V+lx96rcn4oNImCcxM6XPgZnVnvFaHWw98ljkz3cTNRiaVdB0YHf8eKMYCdBGIlV2U02zNhJjmRoJlyfRuER7h14jYQyQ25OIkOGUr0ykDeu6ePT2q9l71yt59Par562BgGhy049EpBFARN4qIp8QkVOjfoCINAH/A3xQKRU2SvmJi75TZxG5RUQ2icimnp48rd2cB7wz4yYZY1g1pNYyW1bpx3TjEqZGIltPYtLHkzDv5Q5cq1kdvzBMDMRTeDMwEqEzNtO3yU3LqTAU3Ui4B9wFos/Vkcmm+IAbZiSM12RmtaaozzQ9TJOowdCsgqb9u6Fzrf49qpGoqgv3JKaGdDwhG0/CHZOoqNKGwltQZzyJhs64kYjgSZRSIz0AHrkNDv260EeRQBQj8XlgVEQuBP4XsB8dkE6JiFSjDcQ3lFLfdzYfMzKS82iurEPAKa6XLwd8fXel1D1KqfVKqfWdnZ1+uxQl3plxs4wyJo2ptcxMayWMJ1G3QA9mmabAGk+ipiUekzAtOeoW6MdqZ7B0xyUmB6Fhib6pM1zlLXDGNtrjYyRWJqfphuAeWBdU6EGzXzXHtxsj4SfzmQHK/N+1znmYOBnps71EzX/POE9+7IT+zhZfpJ+HxSSqPZ5EypiEMRKZpsC6spvAv6DO9G2qbYuf8whGopQa6TF2Ajb9Czz3rUIfSQJRjMS0UkqhYwafVkp9GlIn5DjZSl8CdiilPuH60/3A253f3048vnE/8DbRXAoMGFlqvpA4M1Y0ySgvXrMytasay9zZl94HTrh6Q9W2Zp4CazyJhee6PAlX3yaAKqcJnjsuMTGg9WrTUTVXqFldxd3gIzfNTESe0boH1gUVQ0ypSobcnl19h07HdAfuDUbqMLPamiYtk0xk5klEDYZmHDQdcOIRixwjEeRJTPl4EmFy08Sg3r9+kTaosyHV3H7MTOpzXO02Ej4FdabaWsTlSaSWmwpRfJYxfdv0Yza1PnkgipEYEpE7gD8FfuxkN1VHeN0VzmuuFpEtzs9NwF3AtSLyAnCt8xzgAWAPsAv4AvBX6f0rpUFsZvyPV1EhijUrImiZta16xp5urYR7bYDa1uw9iYXn6HqHyeG4rOLOboJkT6KmNWWWUNqBxfF+bSi8noRZa3s42tzCPeAukCH6VTP11VXxATesoM7MYo3cJBX6XIxn5klEDYZmHDQ96WQ2LXKWwfUzEtMTetD2xiRSyU01LY7BVum3YHEvXWrw8yTGeuKTgjTkppLKROrdqh8zTH7IF1F6N70R+BPgz5VSR0VkBXB3qhcppX6Lf5wB4BU++yvg1gjHMz+I0gHWTcuq9D0Jdx+empbMPQljJDqcbOihA8mehGmnPRXgSQQYiYxW6PK25DCY2WjEFfLM+9/94E4WzA4yLC2JA67bSLSdlvhiM4s1xhH0ucjQkzDHEyUAmlGevPEk2tfoQdbPSLivF0NlbYoGfy65CeJV0VExE5lqj9w0MZAY0B49Hk9USCO7yf0dz1Vr7YyJGYni8iRSGgnHMPwPun4BoBe4L69HVQ5MpGkkWlemv8KbV27K1JOYdHkSoHX/QCPh9SQcI+HnBXU/yqcfHGRsKjGGECuYCzQSnpYcBjOgTEdfRjU24H7n/8DsCla5PzPMk/DKTaC9vQw9ibzTv1uvw1DdECw9TvkZiZDAtVLOQO4xEukQ8yRcRqLZlQZrFiIa69EdkSGt7CaY++IzN76NANuf023R289M3LnPMRKjx3TmXEVx9F9NKTeJyLuA7wH/6WzqAjbm86DKAnOThrUJd5NJrcRUjjyJ2MphZ+vHwf2uwLUrBRbiA/T0hB5calt1gZp3oB05Bt96Ca+a/KbvR4YGFkeDPAnHULmzsaLi7ttkSEduAh28LmYj0XaG/r2mNbonUVUb3G5kagRQ2RmJKT+5yUnZdqczJ3gS0eWmQhKUfjt+/5vh17cl7qyU9iQqa7WUOnKsIMfsR5SYxK3o+MIggFLqBeK1DZZMic3yQ9qEu2lZqW+KdJqoTXpjElnITdWNWgaoqIrLTVIRlwm8cpNbTqvv0Pu7+wX17wIU19c/6fuRoYFFb98mg583ExX3gkOGUCPhIzfVJctNuS7kypiB3dB2uv496FpIV25y719vjESaGU6TPnLT4ot1wHyXk9MyPaH3M4YojeymQuKXfjszNU7dZA90P5o44Rs5quN9y1/qPC+euEQUIzGhlIolv4tIFQH1C5Y0MINoVE/CuNrpBK8nBvVNXlnjeBJZBK5rmp1+/8v1DM+05DC9c7yB69j/1xofbN1N/pz/49yK51hakzjzTxlYDDQSJiaRpiehlL8nUdOijWJUual2QUIKbCFbSicd6/DhRCPh51XGBn3XNRkmN7kXKKpvdxoiZig3uY1EbQucvgF2fksH0t01EpBWdlMh8fOGF1U4gf3xPjjhqtUw8YhTnfVahoonLhHFSDwiIh8B6kXkWuC7wA/ze1hlQLqB60xqJdzpjLWtehBIN0UREls5t6yIxyRMPALis3hz45qZqvEkIDHzxTESguKzl/Wkl60z1qNlHvcs3n0M6XoSUyN6IPQaHZHgoLuv3JToSRRNIdfAHv3YmqknEWQkXEZFKrQc5JabZqbge9fBgYeDj81PbgLd5mT8BOx5IP6eXrlpqrg9CT9veEmF6x7o/m38d5P+uvJ6/VhinsTtQA/wLPBu4AGl1N/l9ajKgbSzmzKolXBXz5rPSXeWbV4TMxKnxmMSbiNRFeBJuI2Ee7Ad2KMra+s7WK8eT6/FwVhvco0E+KfhRsEYL7+V/wKNREDgeno8VldQNIVcpmeT8STSjklEkJtAfyduI9HzNOz/eeJgGPgeTYnbT71W9+La/jWXJ+HITWadi2KSm8b64KfvSKi490u/PbXGeNMChx+N/6F3q/7/Fp6juxsXUYZTFCPxPqXUF5RSb1BKvV4p9QURKZG1GouYKKvSualp1oPy0MHU+xomXW2fa53YRyZxCXeXzuYV+gIe601cqtIbk5hwy00+60UP7NUz25U36MZ4EaukY+/jnfVDfEU1T+A6ZVzArwOsIchITAWkwEIseF00hVxeIxHkScTSUSNmNyUZicTWHE8/8SAA9/zi6eB4zJSP3ARa5jv7LbDnR/GsPuNJSIVT5FdEctPu+2HbVxMGfr+alndc6HTWXf7SRCPRtxU6ztOSbuPSoqqViGIk3u6z7c9yfBzlx+SgnoVWRqlLdGjoTK+9xeRgsieRSVxi0uNJqBk4sSPRk6is1jd2zJPwkZu8nkSXJuVVAAAgAElEQVTbabDqRr396Kbox+PtAOumujEhBTZSXGAsxJOoW+hfIDY9Bkh8zQWItyhxJKeiKeTq3+1kmbXr57WteoD1LjyUrtzkNSr1cSOxcXM3+7b+EoBGGQ2Ox0wOA5IsHQKc8zZdjf30v+vnbu+xwKvTJWEaO44cTdjsbStzXquTBLLqJjj5gs5iUrPQuw0Wnqdf1LSsNDwJEXmziPwQWCUi97t+fgnkYWX7+cvGzd284q4fJ85kJwaiS02GuhTrG3iZ9MQkIHtPwqzbMDGQaCQgsV34hCswb2boZjCemdQN3FpWORqswN4HIh3Kxs3ddB85xH3PTfrPTmuaEjyJSHGBVHKTX8aOWXDIveiNx5MompbSA7u112aO1VwL3gnD5JCT6OCauIR1gfXzJBxp6O4Hd3JexXOAbmQJAfEYkzknPkPRoguh8wI9mJq+TYYCr06XxJHH9eNIimr/oW6dJdh1pX5++Hc6xjc1rD0JcCrOi8eTCKvW+B1wBOgA/tW1fQh4Jp8HNZ/YuLmbezY+wIPNf8lrq/6VZ/tXc8f3n2X9mcdYXhsx/dVQvzC5O2YYk0Px3PiggSHq+xhPotnVANhrJNzrXMc8iVata1c3xg3c4H5AaU+ifiEsfTHs/Qlc/tHQwzBewRMt/ZyYbaF7yKc627OMaqS4QJjcFJQVNj2aPPv1eBLmuApe3du/O96OA+Jp1xMDiYbR7XkaTBdYv1XgvJJpwyL9ntMTDA/0cHqHNuCNEpeFkr4Pb3M/L+e8DR7523jfpthx1Ucupss7UyPQq6/DlC1hhg9pI7D4Ym2Au38bb6ff4fIkDv4yf8ebJmHrSexXSv1KKXWZUuoR189TSqnYqvYi8vu5OdQiYHoCvnk57H8o8kvufnAnS9QhqmSWi6u1tjo2NcPBo0fT9yRSrZTmxU9uytSTiMlNrka9fp6EGaAnBnXqbVVt8rGbNF6T1rvqRjj6RMoc+7sf3Mns1BhNFWOcUPr/SZqdVjclBOcjxQWMJ2HkGDc1zXqQ9Eoz02OJmU2QdSfYvDA7rZMdzGQBgr1Kb5twcK1zPUkSU0NoqciJR5nA8lgPV7ft02+pqmKeBPh8H+5ry4+z/kR7Gd5EhWKSm45uisfURo+G7zvseBJVtbDkRTouYdJfTUeDpi490SgSIxglJpGKuhy8R2lwcicc+T3sezDpT0HB0cP9Y7Q5axWcWRmvIK2aGY5eI2Go70ivgZqf3JRpTMK8T3VjPB7gayRMMd1AYqGgu+rapGSadTJW3QQo3/Pq5nD/GG0V+vhPzLYmbE88hrgnESkuMNqj/xe/+JD5v72dYI3c5MYE8jNcUyIvDB7QhsIErcF1LUQwEiaTyC8uYfY3M3xX1fUtq3uYVcLmqTU0OkbCNx7jvrb8aFqqA9jLrkjcXl1EcpORmjovDPck1KyWkcz6LsuugGNPwbEndJsO873EGlUWh+SUi+Yg5VNYd0JrrPS/kLA5rEndsrZ62ib0zHZNVdxILKgay8yTmB7TA5RfoM/N7IweLKuz9CRMK+eE3jorkrObIFFumhhMNIJeT6KyJn4zLL5IDzB7fwLnvDXwUJa11dM6ZIxES8L2GNVNCd5WpAZvQdlSEB/ApoZ0wZjBLF3qxhjNAnoS3l5B/3LxQS4DaD8rvlOYJ1Ed5EmEGAmDyT4aPc7Z7GCwaTWDQ4tYNLuTrqDGelND4Z4EwI0+y9dU1RdPdtPRx6H1NOg4Pzzdd6xX30tmpciuK+GJf4bdP4RTr4nv517CdcEZye8zxxRHB6lSwVRIuheUJzw4etv1azj8gDYSqysPAIr66iqW1U9lELh2pZJWrwjf10guZqA2wcF0PQm/itiWFXD8qfDAtWnuZ6jviHciHdirs6QqnBm+VMCKa1KuyHXb9WvY8eNvA7BvRhuYpNlpdWNSLUjKuMBYT3KzQEOQJzE9miw3VdbobQXq3+Q3Wfn947/isnriejckxiTcTA4md3A1RsIveO01Eu7+TUcep+X013BtRRXs3s6j77na/6Anh/1rXlJR1aAb4RUDRx7XKa2NS7Xc5Be/gXg8MeZJXK4fZ6dg4XkxA98w9AI/b4cntm7lRae8bG7+hxCiNPh7r4gsCNslh8dTGKbHoefZ1PudNEZid0Jef1hwdMO6Lm44Q99oLRWjrGsb5mOvO58GRjLzJCBaGqw380REf166noRfHrsp7EsyEm5PYiA+YzXH7pabWj2tt9tO01WmIRXhGy5YxAfafszm2Qt4fmalf7aQR26KRBRPIorcBDp4nUW78Gzwm6ycLns5qhYlfhdBnsTwYT3QuQmVmwb9jcSRx7QsuvTF+u9+izYZ3EuXpkOxxCSGDukZ/9JLdXHo9HjwPWbSWo2RqG+PNc18cqQrlqp9dFZPBn/55ObC9fpyESUmsQR4QkS+IyI3OCvOufnTPBzX3LLtq/BfF6fWko3cND2WoD2mCo6e1hQP+t33xlY2rF2WLMdEIazhnJdY4zTXDZhJk79Y2wSP3ATpexITAzoAPLA3HrQ2NC7ThjesgeHz36Vhopt1r7szuDrbkwIbiaAKboifvyRPwkdugqwWHsoWv8nKmVX72T7p8Tr9YhIzUzrH33RgNaQjN9U06/13O117lr5YXwPTo8HGP1V2UxDFkt1k4hFLXxw3sEFpsF4jAbFU2M89Uxsz8EOqgVFVS7vqLYq1uFMaCaXU/4deS+JL6CK6F0Tkn0TkdOfvW/N6hHPBcLd2+cKqmZXScpPRdl2SU8rg6PiJeOCwd6u+adRM9A6wsTf1qVwOwq+BYCZN/vzkppXXwYpXJC/EUxXmSTjHPrBHn48Wr5FwbrCgYJ1S8MTdeuZ12k3Bx2s8iagt1ZUKL85LR26C/HoSk8Oh/5d3slLFNKdXHqK7yqNrV9ZoD8E9YRg5Cqj4Wg6xfY2RCJKbXNeXiPYmhg/pc9NxniumE2C4U2U3BVEsgesjj+vz2blWB9kh3EhIpfY4DGf/CXRdye/6F7t2FI7OLGRJxYmiWIs7UnaTs2rcUednGlgAfE9EPp7HY5s7zMwvrKnWyBF9QZ/2Kv38ZDx4nbJoavyELmZqWqaNxITPAB6FTDyJmjx4Eh3nwRt+EU99NCRkN/l4EgDHnMpqr4FJldFx4CHo2QLrP+RfeBU7hiZAJWvofTvgwXfqTB83U8M6OJ9uTCJIbsqHJzExAL/+MPz7QngoeFVf72RlVWU3NTLN2edflryzt3/TsNHLPZ5ErE9SBE8C4pLT4ot1AZz5u9+CV+5Fi9KlWOSmI49pA1FVCw3O4D8SkAY7dEgbiArXhPKUq+BNv6GjLXHCeGx2IYsq+4piLe4oMYn3i8iTwMeBR4HzlVLvAS4G/ijPxzc3mJs6LOXMBK1PvUYXv3iC197y+wQZZOKkzsFfeJ7u9phucz9D3QJAoqXB+slNGXkSPv3+gzAxCaWCjcSRP+hHr9xkjESQoX7ibn2DnR2c/aSPwfSQ8sxc9/4Etn4Z+vckbg9awMgQ6EkEyE1+nkTvNviPpemvUT47A0//B3xpNTzxcVh4rn6+/eu+u3snK5e16HO5/mKf4Kd3wmCkkOYguSnAk/DGE4yRWPpi/ejODvMyPa496kw8iWLIbpqdhmNPxv/XKJ5Ek38ChdfAH5ttZ2nFiaJYiztKdlMH8Dql1H73RqXUrIi8Kj+HNcdMRDESTjyi/RwddPUYiVDGTujBo2kZPP35uFFK10hUVOn3yVRuqm2Fk8+n95lBDdj8qG7Ust3kkL6BvIFr0EVzkBy4blgMiP93cPxp2P8zuPKf4sV5gcfgXufa5R2YAXH0aOKykeZcBsUkgga56SBPYkFyCuyRx/Ts8sDDcP47w4/fzebPwq/+GrpeAi//ic7D/+418PO/hEUXxdccd5GQyfXb38EfKhPTX2PH6VlTwpt5YwiLSUz5eBLGI/MaCb/gtZ+XGpWqBn2NFXKZTyMdL71UP69p0ddEUK3EcDcs8B/0vanaI9WLWVb9GKesXZaPI0+LKDGJ/+01EK6/pbnocpFiAtZhTbVO7tSDYHOXzl2OaiTUbNyT6DhPz0B7ntZ/S1duAq3tZyM3ZRqTiCIJmFm8mUn5eRI9m51mc56EucpqPQv1m4Vt+Zx+7wv/MvoxeD0JMyB6b+CxFJ5EZa0ehNyDnFL+FdfgrCkxkBioNctwGgMZlZ3f0e003vhIXL555X/r7+KHb0gdoO/dCgvO9DesSXJTt5aWvFXnlQFy0/SElumC5KYlxkiYxpJhRiJDuQkKKzm5g9agYzKNS4KrrkM8CUhUI97yisuonJ1IXKirQOSi4rr0ieRJ7NQ3nFToFgf9u6IFRyeHtKEwRgJ0Uy9I35MAp8lfFnJTLmISQZhB0wz0biNoajymx5O9CENQi+Tjm3UWiNew+BG08JCRgLxGKLbKXYAnIZKcxmniHUFyE3DDx38Uq74/sM/xQtMxEiPHtAdyxobEnPumpdpQnNwJD90a/h69z+oCLz+8ctOQ01PIm7xYFSA3+U1CQFdHX/bRuGwV5kmkI2V6iS08VEDJ6chjenLhvp4bl/pPdCaH9fn2ynlBxArqCl91nVcjISJfFpHjIrLVte2jItItIlucn5tcf7tDRHaJyE4RuT6fx5ZAlMD1iefirmLbGXoQCgpQJby3MxOoi+dEZ2UkovZvmhzUM88qV9eU2lYtGwQtbO/7Po6R8Js1ezEDtJmtu7O3qmrjg4E3HmFoWuZ/U/ilzAZR45abXJgB0fudpYpJgDa0CUbCZ8Ehh6eO6YnDyGBvrDX5kYNOPKv32eCOql72/BhQcPprAE/bl28q9i57C+z4RvBEZXJYn7fOECMx6fEkvEFriHsSXrnJyG/ea3jRWrj87+PGJmYkfDzYySzkpmJY5/r4Fi37uQ1r41J/uckv/TWMWCKHS90Y7clsZcksybcn8VXgBp/tn1RKrXV+HgAQkXOANwHnOq/5dxGp9HltblEqtScxNaYlA6PtLlitH6NITjEjsUDfDK2r4q/LyEikITe5++q4Py8dyWlqWBuIighfRXWIJwHxgdib/mpoWpZsqCcG/VNmA4/BMVReKSZmJHw8icqacMnD60nE1rdONpzfflZ/bktF/POXVBxnRDXoeE1PxAbKu+/X9SidF/quifHtXXU66Bv0XZrlMM0aBV6SAteH/AewoJhEkCfhJUpMIhtPolBGYtZZU8XrqQXJTWkbCY8n0fcc3HMKPPuFzI43C/JqJJRSvwaiimo3A99SSk0opfYCu4BL8nZwhukxra1W1upZpp+l7t8FKGh3eRKx7SkYc3kSkHjTZhST6IDxiEbCm3mSSZM/91oSqUiKSXjqQIyR8Ka/GhqXaZnFnaZqlms1a3ynPIYgTyJIbnJqJPzaKBiS5CZnYPKRm/YP6yaBraIHQGGWpRW9/HziRXqHKJLT1JgO1J/+GhDxraTunXLOdZBmbToIBHkSNa36HM1O64nScLe/FBLUliOykQiJSUR9Dz9iRqJAclP/bm04OzxGuHGpVia85ytdIxGrG+rW38/D79Of1/1o+OvyQKFiEu8VkWccOcoIzV2Au5rtkLMtvxipqf0sPTNzr9FrMJlNRm5qOVVLOVGMhPFSjJEwF1VVnZ7Bpkt9h74AU2mxkz4V3Zk0+XOvSpeKsJgExI1EoNy0FFCJ34G3rXgqggLXYZ5EmNQEehCbiiY31TTq2EuLaCO1qOIENTLNC9VrddzD1ImEceAX2hA5UpNfQVW/cgbWoPhU31b9fQSdt1hrjkF9DmYm+eRjI8lLvAa15Yg6wPsF/g1ZeRJmPfMCeRJ9joLuzTCLVV17vIkhk2IccUirqtXX5chh2HWfviaqGnSt0BxTCCPxeeB0YC16USOzoJHfVM5XcBWRW0Rkk4hs6ukJX4MgJWYQX+h82X5xCdOzyaROVlRBy8qEgrpAxj2ehLmoMpGaILHJXxhhnkRaRiKNtglh2U0Qr7oODFz7FNQZIxFZbsogJhEUtDak4Um85WVrAWh15KauCn19XrnuIr1+QBRPYvf9+jOd5m5+BVUnZ53vNsiT6H1WX2tBhYeu1hy/3LQZgOeGW5OXeA2UmyKu0e4X+I+9RzYpsAWWm7xrQBgaAwrqhg/pzDdvAWoYTcvg5C745V9rWWvde/WENWpcK0fMuZFQSh1TSs0opWaBLxCXlA4BrhVtWA74BgmUUvcopdYrpdZ3dmbQQdKN8STMDN8vLnFip+737v6C2yKmwbpjEhCXmzI1ElGrrv2qYTONSaTrSQwHGQkTkzgVX/yqrgf26s/3W1rUDz9PQiktN0lFbNYcYzyiJxExJnH9xTputaJhAgHOa9Yy12XrLtZG4sSO8NRVNat7H628MeZp+rV9Ga90BvlAI7E1OLMJEiYMD/xeGy7TWA5cizlVVAGSudwETuDf55qLeRIZyE3VBZaberfpyY530A/q35Qi/dWXpi7tQQwdgFd8Dhav1/Jg3/bMjzsD5txIiIi7zeRrAZP5dD/wJhGpFZFV6H5Rf8j7ARmt2ngSfkbi5M6kIpjd00sZObaTVbf/yH+tZcPYCT3rMW57+xrdvyVrI5EiDXZyKFnuycSTSMdIuD2JqvrkRXwueDdc+4XEjCs3flXXJrMpLGbgprJaV8S7PQnTK8tILyOuFtNRPIk0spuobgKp5NbLO9h71yv5x5c730HLCn2Tq1md0hvE0Sd0C+wzXhPb5Nf25b03OXMrPyMxelz/ePVyN6524TVj+ny7jQQ4MpdIfAlTNzEjEeE6rm0JjklIRfD1EEZVgbOb+rb6n9+gqutMjITxrM96s25Fvkh7qRyfW8kpr6WKIvLfwFVAh4gcAv4euEpE1qKlpH3AuwGUUttE5DvAdnR/qFuVUvnP93LHJKQiuaBOKe3infO22KaNm7vZ9nwtf9cwSrsM0N0vyWstG0whnaGqTmdHpbu+tSHIk+jboQdBc8NNDibP0DLxJCaHE1IjvYvaJCwkY4zERL9TQe1h4dn6J4iGRc534DISg3uD5akgapoSjYQxiu1n6YDj6FG9DOvMlD7WqJ6EWScgRG5CJLF/0+B+J7OtWXsSoA3B8pf4f9bu+/UkYuWNCZuT1sSYmYLf4j9ZMEHriJ7E6oZBplUFPbOJdSgxmauyNgeeREBMorop+gTATSHlpukJ3bng9JuT/1bfqa/hJLmpOzjTLIiOc/W19NK79fO20/U9Zopx54i8Ggml1Jt9Nn8pZP87gTvzd0Q+mJhEfYce2LyexMhRfYG7Whvc/eBOVk/pQXBl5WH6ptti7nmSkRg/kVzFevVn41pvuvh1gp0cgq+vg0tuh8s/Gt/mvYEz9SQczThsBb4N67oSB81MjGBFVWLVtVIwsA9OCViwJoiqxkRJx3iLC84CfhyXw0wPrFSL3tQ0a09kelzLHCFyE5DYv2loPzQ78lrjYi1bhsUldv9QG5D69uB9QHtMNc3+nkQsqBrBSEwO8PKuSXoPL2CWuKSV0MW4stbfk6iqi9YSo6Y5ealUyLy5HxS2mO7k81r28fMkKiqTOwfMTutxJGohneGiD8D5f+FaE6YCOi6Y8+C1rbg2M77aVv88fRO0dslNh/vH2O+sjLay8kjCdkgsfNr8/B56pz2DyanXBM8kU2Ga/LmNRN92fRObPv5K+ctNpkV0up6EIzeFrcCn3782HijNVE5rdBXUjfVpIxU1s8kQ5klA/AaOVVtH8CQgnuEUJjeBx5M4oKUmw+L1wRlOM1P6u/Su5xxEXbu/kejdphMcvKvMJRxjfMJwau1JatpWBHcx9pObpnwSI4KoafbvAjsZYenSIApZTGdqUILkvIYliUZi5KiWGdOVm6Qi2YguWqs9iait8HOAXb50/KQe0Coq9QA1dCDx77HGfnEjsaytnkP9i5hRFZxaeThhu3e2XT87wJbjjQxv7g5fQjMqpsmfuxOsCWQdf0pfkKZdtt9NnG5rDldMImwFPsDRrxv0azKV05qWxZvNDaaZ/mrwLmEaMxJnAhKXAlK15DC4C8IaFoXLTZDoSQzu1+2gDUtepFMax08mtxkZOqA9FrP2SCrqFvp3BB45rA1TaO2Hy6scOkT74rN59F0BHpuf3DQxGN0LCMpuSqcGx0sWclOoZBqF3q1aEgxo1keTp+o63RqJMDov1E1CB/dHrx3KEutJTLhu1qZlyTGJ3q16kHS5irddv4aq6joOzS6KeRLGPffOtltlmL6ZptyuMOWtuu7dFv9934PhenE6Tf5mpvQM0rmRU63AB8TjEpl6Em5vLt0aCfcxJHgSzoBd3+Hknjs3cJSWHJBcNTwV0ZOYGNDnutnlScTiEj7ehGljHjUGU9ceL9Z0M3LMPybkpqpWD/4TA8EtOWL7BgSuIxuJFv9W4VNZyE0VVfonzewmv+r1WLpvVHq36rhiUEdis9a1od9Z2z1ducmPAgSvrZEYP6nbO4O29GO9ib2Njjymb2xXvrnJNjkqy1lddYDlbbUx99w7226rGKZfNed2hak6T/+mvm3QeYHO0d770/hg5lfRnY4nYQZax5NIuQIfxGfXmVSTg3ODHY8vcwq6JiUdqgPkpppWPcvzyk2pYhLeJUynR3UGVZAeX+e0CzfdX90pv4sv1o9+ktNABkbCN7vpWLjUZKht1fn7k4Phs1y/mMTUUPSJgDvw7yYbuQm015qmJxEkmX7w21vCsxTd9G0LzxxrXOJ0DnA+Z+tXtBEOe01UOs7XY9EcxiWskRj3eBIQnwVMOW29Tb94FxvWdfHil76Os6v28dtTPsKGpXrgcc+qa5mkXiYYUE25XWGq3tMJtm+7zpxYeYNu52CC8X5yUzqehKdLZ8oV+MDlSWQhN4Ee6Ab26oEwXYMTJDfVtibqxaZNuClQDMLrSQQtOGSodeQmPyNRt0DLSceeTH7dwB4dN2qKuIZAfXuy3KSUYyRSeBKgz4eRKsNmuZW1/l1g05Gb1GzygJ7p0qWGDFanC5usRfIqpka1ZxCWqdS4VMuGY7066/DAL2Dte3Kz7kV1g+5GPYcZTtZIeOUmiAdOjz+lMxN8jASgs4lu+i+9NvY3LoFfvIcPX7syNts2Td5GpSW3K0y5O8FODmktu+NcWHWjnlkeeFj/ze8mrmmJvrymT5vw0BX4IJ7xk7En4foOBvelLzWBjyfRrzXk6kbHk3DFJGrbkus5vPjJTUFSE+j3nJmMx7PcgWvQs0G/gqiBPdpritJMEZyYxAk9ABsmB/VnN0YwEjWt8WMMNRJ1yZ2D0zUS5tgS3iMLuQn0d5BmdlOqyVpCIoYfJ3YAKoUn4WrNsflz2sie/660jjOUzgut3DSnTPTH5SZvW4gjj+lHs6iIFxHdP/8dz8HaW+Hp/+A1db+JzbYXiB5UbnrRubkJWhvqF+pKYaX0TAX0inkrrtGu6M5v621+N2DnWj1AeZfx9COTRWFcnkRCe+uorrzbUKfTItx7DF5PorbVWRTGMRJKOYV0KeIR4ONJjIa3TjeTjp6n9QDhlX7az4b+FxIrv0F/J1GD1qC9LDWbGBQ2hYJRPQkjI4XJTVV+KbDpBK4Dmvyl0/LFj+r05SY/ydRLqDQca8eRvCpgDNOa4+RO2H4vnPWm1JJmOnSu1RMos1hanrFGwk9uMsHrI4/pmV2qWVltK7z8U3pA6N0Wm23/7D06T/2Ss9K48aNgmvxNj7rS8c7V8sPSS+OuqJ9mfMG79Ex1y7+n/hxTa5COJOAYic3HZjILEMa+g0P6Rojas8mNNwV2ciCebdW4RLfsHutzmvtFuHmTUmBTyU1t+rHnaSfLyHObLTxHe6juti5KwcDu9AoHTf2NW3ocTdNIGFLGJLKUm8xrDGpWf0cprq3QiUYGcpNbMg0i1Nvo3aYlwQVnBO9jPIk/3KX/x3XvS+sYU2KC170R285nSXkbiekJfZEZI1HfoQOSxpM4/Fiw1OSlolJrhcZ9h3hQMVVhVLq4W3P0btM3sRlcVrkqdf1u4qZlcMbrYOuXUrvqmaxB7Myw79s2FF5TEYSpWD32lJ5pZ+pJzEzq4Dc43qIzcJsbePRovE14KtKVm8z11LcjMbMJPei94359Xv/ui9+JD3omGyodI2EKK93B65iRiBC4NnGjuoXhrTEqPdlNsQE+CyMxPYZO0w6+tlJmIlXVZ9S7yUziPvXGtakTMbz0bdWeYFh8wXgSxzfD0sviyQq5ovNC5/3nRnIqbyPhBHj/5ZFjeqbyz79ktGaxTsEc6taz2WURjQToYq2TbiPhaROeK9ydYE9s159rdOyVrjWegmZ6696nB84d3wj/nCw8iQMj/jp/yiyvikp9kx12+uZHNBLuGednfu0YeeNNTLg9CdOn/4g+f1FkALPCnztwHSY3GYOkZhKC1mbQe6xff2b7xK74oJduZhPEryu3kTByU5SYhDknqVIzvXUS5rqIGnfyMxIR2nqkLN70ZjdNj8P+X0Q7JuJexasXPMeHG7/qn4jhpXdruNQE2ss0Xvy690Y+nsg0LtGTAGsk8s8vNusB/cBIbWymsnO4meNH9sJRs8h5mkaif088yOdtE54r3P2berclXrSLL3JmkRLclrjrCj0b2fK58MrNTPr9OzJMXaP/etSRsrwal+nWBxAp/dU74zw8qg3mT7c4cs7EQHzW7O7SGdWTgMT+Q9Oj4XKTu0jOZSTMoDdGHQdnFrO68mB80DNGImhBJt/PCZKbJNr/ZYxEqiIvb0wi3cWCqn0C1+5OBwGkLN70yk3bvw7fuzZaC3+HDeu6+Oya3/Cehu/x6AfWhhuIiUGdpBIllbVxqR7Mz3x95GOJjIiOS8xRhlNZG4nvPqr7Dg3MxgfBIzMLGD1xQEtNlTX6y4iKWbhowCmeGT/hX1qfLWYAGNwXz2wySIWWnBo6gytuRbQ30fMMdP8m+HMykZscw/THl5+TvitvcKeARqgq9c44R5WWTu79lVZSVP8AABvqSURBVBOvmeiHOiM3maDi81qSihKTgMSq4SjZTQaX3OQe9F6YPoXVVQfi29MtpIMAuem4I5tGSLeM7EnU+RuJqHJTrU/geuhgys9OWbzpzW4yGWNRl4gFPUkyXqtfWrIbM3Fx9XEL5Mo74fovZ7awWBSWrNf38Ry05yhrIzE1omdgAyo+CB6bWUi76tVB60UXBVdV+mEuHhOXGD+hM6eCFn7JFDM4HHIG+HbPwicv+1f4o5+Fv8dZb9Yz3s2fC94nNhiksVCKI8NcfeGZqWsqgjBGonFppDbS3hnnqNKDyNCwk/3h9iRqmrRn1Ot0So3qSbiNRFS5CRI8Cfeg98LMCk6r7KaCGb19YI/2ANMxyCYrzxuTiCI1QfycpPIkkuSmiAsOxT7HR26KGYlTkvd3SFm86c1uMoO4yUCKwsnn4+nkqVYNNAktIccc48w/SowP5por74S3bsqsg26alHXvplMb9ezIbSSOz7bTUjGq5aYL35PeGy5wVq6LGYmTuQ9aQ7zJ36FH9HPvEor1C1Mv0lPdAOe9E578pO6V5DejmxyO3unT0NCpg/91C9iwrjmz1F8jCUWMRyxrq6fbZShGlTbsK5pVPEU0IZNnadxIRE1NdC9hmkpuqqx2ajWGE4zEbdevifX12jWznFqZYnVtL++5/mLYlWZmk/kcbyfYkYjV1uCSmyLEJGan9LmUigzkJmeSkWAkDgESaqDMtRPYZ8krNxkjYTL+otD9W+cYmyIYCSfWFbXYcZ5Q1kZiwzmN8Hyi3HRCnEFjZjK9eAToWWDT8kRPItfxCIg3+Rs6mJjZlC5r/wqe/IRuGHalT4f2TCpiz/0zfd6ykdhMvUpEI+EefCHuSbz14oVO91GVOLtvWBKX2dLxJMza26nkJtCfNzWSYHzdg96uYS1D/cMVwqXrumDTHlh2ebRjceNt8jd6LLiux0uqlQINsXWuJxMD+FG/Y6nQ15E7JjF0UHs8KeSYpHU0Eo7Lld00Mxlv49K3LbmJ33VnsuH8Bclecfdv9Tk89dq47BTEcLf+X6KkF88jylpuurBD63lNrR0xSeTVV6yP75CukQAtOXnlpnxgbnB3ZlO6tK7SC6c8/Z/+6bCZNGCrbtDB82wwM7WIPZu87UKamvQM+fIV9fF1DNyehPFUIHpMojoNuQm0EW9amjQImvTL+/7/vwDg0janT9XQgfSC1rHPaU+Wm6IOYstfAq/+Lqx4efh+Zu0TIzkZjyqdJo7e1emGDkaTbcKoatDHpJQ2EGoGmk9h9sTz/O/vP5WQOrvlx//M1L8vS16k6fCj0HWl7s82dDBx1UIvw4d1TCvT+61EKWtPgomTUN3Erz90XXxb3w7Ygr7RUs2w/Gg/S1dZKqXf30hQucakwaZKx0vFxR/Urat3fEMX2rmZzMCTyAVGgkjDQ0qYcQ7shS+ij3/Cx0g0uY1EmjGJmSktvYTJTaAHE28rcDdm/ZK+7U6L8FlozaDo0t0JdmpE/0Q1ElIRLfsm5kmYrD2TmZSGkfCuTjd0MHyVwijE2oWPxzOazngtFZs/w5LZAwyyMrbr5ZWbqJ4ehOe+Betu1RtHjunXnX+LDgSDDl6fdpP/52WyBOk8oKw9Cd0Bti1xm5nFLr00s6BQ+1n6Zhg5mj+5CeKDmzcekS5dL4FF6+CpTyVnSmTbgC1TOi+A674EZ70xs9ebY54aibcJd3/PxpOorEl/TYRpV/plGNf8B1z/lfB92s9JbJGSiSdR75KbjBwWNXAdFeNJGCPRt12fz6heGCTGdJTKjSdRbYzEaDwesfq1AKyp2u/aUXFRtdO+ZpvrOzHyUtcV+h5AwuMS1kiUIX4Lv9S0wIpXwJoMByiT4dS3TfdWybeR8GY2pYsIXPRBfeN7C5Ems1gUJttjOv/P08uqcmNeNxXgSZg02PqO6BOBmmb9fqZAL5Xc1HZa6kF/4Tm6YZxJmc4ktuSWm2J9myIGrqPilZt6tujWEOlMotzZYRMD+lzmQm4CbbhPPq+966WXMa0qWF0ZXzxsZeVhOioGeEGdpj0FswZ492+1l7ToIn187Wf5r/NhGDkcj5eVEeVtJCZ8jIQIvOEXcLbf8twRMEbiyOOACpccssFkL2XrSYA2iA2LYfOnE7cXypPIlqp6QBxPwrWWhMF4EunOhCHeXjyV3BSFhWfrYzz4q/RahLsxRkLN+vZtyqjJohe33DQ7owdZ0xoiKjXN8cB1hPTXaMflWp3u5PNa2q2qZazxNM6uPhjbbX2V9iKOXXSXzrzbdq/+Q/ejsOSSeJr7kpClZafG9KSy2XoS5YV7waFc0bRMD6xHfq+f58uT6HqJ/sk0s8lNVa1O993zYzjxfHx7NstLFhKR+Op0Rm6q85GbosYjIG4kjKSTSm6KwkLHC9z3U93IMJN6mvqF2kBMDCYZiZyswgYuuWlcNyWcHk2vyBS0h248iVwbiSlHbmrX8b/mrgu5rPVILJHhJU0vMFXVwpUv/yM47VWw4+t68nD8KR20Nixeryvxhw8nf5ZZLdF6EmWGn9yULSLamzicZyNxxmvgTb/OXabFhX+pZ7NPubyJNFcOy8msNVeYduG+noRLbopKzEg4A3EquSkKRiqcGEivRbgbd/8mY8AcuSll76OoxOSmiXi/oEXpGonm3BsJ482N9ep4Qdtq/XzhuTSN7+fRD13G3rteyc2de6lefrk2wuf+mT5Pv/8H3Yl32RXx91vsBK/9JKdcrlNdYuTVSIjIl0XkuIhsdW1rF5Gfi8gLzuMCZ7uIyGdEZJeIPCMiWeZRRsBPbsoF7Wflr29TvmhcDGe/FbZ9mZ88tpkr7nqY8dFBvvFUX6TBPmezVnJkbEy78IkBLZe4K+frF2qDmI52X50HT6KhI26oMvUI3UZi5JgOKDv/a8reR1Fxy009T2vJZmGasTCvkZDKxFTkTDDfgSmMNJmEHedp7+rEc3oi2LddB6fBaVmzCJ76DCCw7LL4+y1a63Qg9jMSjidh5aac81XgBs+224GHlFKrgYec5wA3Aqudn1uAz+f1yGam9CCSjzoGd2+XfMUk8sEld6CmJzn80J0c7R+mTiY5Nl4dabDP1aw1Z8amutFJge1PXkpVKuCmb8K690d/P+NJmOBwLmISEB9sM8lsgngq9Hhf0trWKXsfRcUtN/Vs0bGUdHsS1TRrIzMzqY1E07LsvWBjJMx6zzEj4cTp+rbFPXpTqFhZrSdDakYbE/f9Wd2gU8odI+GerHz2/l/pfazclFuUUr8GvCu13ww4kSPuBTa4tn9NaR4D2kQky6lGCDGtOt9GokQ8CYAFZ/Cz2Zfzptof01WhZ8wjqj7SYJ+rWWvOJJJqlyfh12n0zD8KXzjGSz7kJogbiWw9ibETSYV0KXsfRcWdAnt8S/rxCEhcnS4X6a8Q/w6OO91QzffZtlp7O33b4PDvtNey5JL46859u37scklNhsXr4egmNj51KGGyUj9xlFFVy8btw8mvmecUIiaxWCl1BMB5NFOfLuCga79Dzrb8EFvrofQ8iXxq/3f3v556Jnh/47cAGHE6qqYa7HM1a82ZRGIC15MDybUwmZCPwDXE4xKZxiTqPXKTq0bCW4meVpNFN5WO3DR0UAd2041HQGKTv5wZCec7OLFdt8Mxqc+V1dC+Rjf6O/yozsRyJ2B0XgDXfgHW/23yey5ZD2M9fO1nv0mYrCyuOMGxmXbu/tnzya+Z5xRTxbVf0rVvH1wRuQUtSbFixQq/XVIzYapGczCAeGk7Q0saVQ05bxVs5BhzARs5BsjJOtpjzWfywMQVvK72YSDeBynVYO/tnwSZzVq9zfrc29OiukkXNDIbumZBZLyeRK7kpnPfpmMImVbOx2ISRm66JuHPob2PomLiOUecNVYy8iRca0oMH4IzNoTvH+m4nGtiZjKW2RRj4bnaixjrg/PfmfzaC/7C/z2d4PXi0WeBuKexuLKPY7MLOTyQ5mRlHlAIT+KYkZGcR2dqxiHAPb1YDvjkooFS6h6l1Hql1PrOzgwXGI+1FsiDJ1FVq/si5cGLyJkcE8Bt16/hC5N/QoVo+zyi6iIN9rmateZMIomlwAbITenijUnkSm6qbYULbsm85XNFlZZyRo5qCTUfzeeMJxEzEmnWSED8/A3s1UV5ufAk3IbaZDYZOs7THsv0aHqNEzsvgIoqLm3an7B5ccUJjs0uTH+yMg8ohCdxP/B24C7n8Qeu7e8VkW8BLwYGjCyVF/IpN4Gu4hzJ/eHnTI4JQA/qr+HXv/gvXlrxOxoaW/nYddEG+1zMWlO2h46KSYGdrsyNt2ikjLEcy025oK493lQy19XWEI9JDOzRg3sm7e9NdphZGCiXchMk90hze2bLfGIPge9ZBwvP5YaZbu4aqnQmZIolFX300pH+ZGUekFcjISL/DVwFdIjIIeDv0cbhOyLyTuAA8AZn9weAm4BdwCjwjnweW0xuypeRuPYemJ3M+dvmTI4JYcO6Llj+H/DgO/jsa98SryuYI3IikZjAtZpNzm7KBNPu2qzWF2ExpDmjfqFu7wH58STc6cOZSE0QbwaYSyNRWYtWqVWykTBLjDYth5Y0P2vRRSze8yM+9trzuPtnzzM6cIxameKS88/n/BxIuqVGXo2EUiqot8UrfPZVwK35PJ4EAuSmpD70mcxiIbHCN4fkSvtPSef5euWrUsXITZAbuQni/ZuqGuZkRbDI1LXHl97MdXM/4AfP9HKz8/tXXmhlwebu9O+Jmjx4EiKxNSV+frSVj373Ydf6EWewoareP4MpFYsvhm1fYcMZig0XXa3bkHwNzl+TZdfaEqWYAtdzy0S/vsBcs6R8B4VzQc7kmPmOO5slV8kJNc1aQiwmqQkS06xz7Els3NzNHRt38MrWCqpklseHT+GRTO4JIzed2OEUMmYYS/RSVc/szAR/89MBhqa04e7uH+OO+7az6GX3cPn6y1K8gQ9mPZTjT2kvpIyrraGcjYRPS46woHAxDcI5kWPmO+4Osrn0JCB3mU25whTUQc49CXNPTFJNFRNsnz6NsdkM7glz7qbHdE1IrtZ9r6rn0FhdzEAYxqZmuG3TMh69JoPU4s4LncrrJ+GMmwuybGnOFI0cUL5GYiK5uV++g8KWOcTdcypDI+G9UTd21tAJuctsyhUmkFzVkHl79QDMtT+haphRFRycXZywPTKV1TqOk6vMJkNtCy/0+f/PGd+31Q3Qfrb2JCDuScxRtXWxKRrl2+DPx5PIWRsDS+FJ8CTSl5v82oM82zOr/1isclMe4hHm2p9U1eyYWYVyhoyM7gkjOeXSSFz3Rb5Y8Ve+f8rqvl18UTzOM3JY99hyB/DzSL7T3NOlvI2EZ/DIWY6+pfBkKTf53aiDM05GU7HKTXnIbDL3xC8n13P/+MuALO6JmjwYiWWX8cbrr839fbv4Yl17MnwEhrrnVGoqNkWjvOWmzvMTNhVzULiYNMqSIEu5ye+GHHGqz4vWk8iDkYjfEx/m8PAYXdlce/kwEuTpvl3kBK+PPTnny5bORZp7OpSvkVj9Oui4IGlzMQaFi02jLAmylJv8btRh5XgQxRaTyKPcBDm8J/JkJCAP9+2itYDouMTIYWcN7LlhztLcI1K+ctNVn4Dz/qzQRxGJYtMoSwJ3CqzpQJoGftLjpDjGodjkJrOUbT6qrXPExs3d/O7gFABv/+6Rwi5IFYWaZl2gd+Rx3YplDj2JnDVmzBHl60mUEMWmUZYExpOobspo3QI/CeOqs0/X/QCKTW5qXKKPyVt1XCQYT/jjdTVQB1sGmvlDKXjCiy+GXfcBas4XGyomRcMaiRKg2DTKksDEJLIopEu6UZ/Z4xiJIvMkalvhnbuL1pMwnvBQbQOjqpYB1QRFWH+UxOKL4Llv6t/LcLEhQ/nKTSWEzbpKn43b9KJSO/srcrfmRrEW0wE0Lc3deuc5xni8Xxt7FR8eej9mVYCi94QXXxz/PQ25qajWes8B1pMoAYo566oY2bi5mzvu286NrdUMqcbcBfqNkSg2uWkOySTLznjCz82s4rmZVQnbixp3sDpiCux8TDKxRqJEKCaNstgx8saIqmdwVscmctJeJWYkitCTyIB0B/xMB8Biy9aJTG2rXkBscH9oryn3eawQYUYlrpVWjK190sEaiTJmvtZeGBljSDXQr5qTtmdMMctNaZLJgJ9pb7OS9oSXXab7OAX0mvKeR6+BMBS9tBaCNRJlynx0iw1G3vibwb+hT7UmbM8KU5SX4/5IhSCTAT+bLLuS9YSv+pRechX/SZXfefSj6KW1EGzgukyZz7UXJtD/5PQ57JvRA1Mu5I2Nu6v5p+m/5YJvNpV8QDKTAb8se5vVt0PrSt9eXuZ5yrcoBWktBOtJlBHumZC/U1zabrEhH/KGDoZvZWzqKgAGS9zzyiStumRjCzkgaFJV6RODAKgUYVap0pLWArBGokzwyktBzJdZYa7ljVJZayQqmQz4JR1byJKgydOMUtRXVyadx0JWSOcaayTKhCjaabnMCjNhvlW9Zzrgl2xsIUuCPK8uV2xivhpOayTKhLDBTGBeXty5ZD5WvZfrgJ8JYZ7XfD+P1kiUCWEzoUdvv7oAR1RalLMebylvqc0aiTLBDnLZUc6DRLGSjzqfsPec7x5DEAUzEiKyDxgCZoBppdR6EWkHvg2sBPYBf6yUOlmoYyxFgi5yO8hlT7EOEvO1KDKMfNT5zOfaoWwQFVAhmPcP1kZivVKq17Xt48AJpdRdInI7sEAp9eGw91m/fr3atGlTfg+2RPDLYJpvmRbziVwM7uX6nV9x18M5l0/z8Z7Fiog8qZRaH2XfYiumuxm41/n9XmBDAY+l5JjPBXK5pBi6dAYVZ6V7LOX6necj22y+ZbDlikIaCQX8TESeFJFbnG2LlVJHAJzH4myQX6TYizw1uRqcsyVXg3u5fuf5qP4uy4ryCBTSSFyhlLoIuBG4VUReGvWFInKLiGwSkU09PT35O8ISw17kqSmWmXeuBvdy/c7zscaKXbfFn4IZCaXUYefxOHAfcAlwTESWAjiPxwNee49Sar1San1nZ3AL33LDXuSpKZaZd64G93L9zvOxDnSxrS1dLBQku0lEGoEKpdSQ8/t1wD8C9wNvB+5yHn9QiOMrVWwGU2qKpSguVynJpfqd5yJon49ss2LNYCskBcluEpHT0N4DaEP1TaXUnSKyEPgOsAI4ALxBKXUi7L1sdpMlHYopG6gcU1ehuL6DciWd7KaCpcDmCmskLOlSroNzsVBOqabFSjpGoiwrru0gUd5YSaGwFEtcyBKNYquTyDvFkgJpsZQr5ZqRVaqUnZEolhRIi6VcKdeMrFKl7OQm6+r6YyU4y1xRqhlZ5UrZGYliSYEsJmxjM8tcY+NCpUPZyU3W1U3GSnAWiyWIsvMkrKubjJXgLBZLEGVnJGD+ubrZxhOsBGexWIIoO7lpvpGLlF4rwVksliCskShxchFPsI3NLBZLEGUpN80nchVPmG8SnMViyQ3WkyhxbPWqxWLJJ9ZIlDg2nmCxWPKJlZtKHJvSa7FY8ok1EvMAG0+wWCz5wspNFovFYgnEGgmLxWKxBGKNhMVisVgCsUbCYrFYLIFYI2GxWCyWQEQpVehjyAoR6QH2Z/jyDqA3h4czH7DnJBl7TpKx5ySZUjonpyqlOqPsWPJGIhtEZJNSan2hj6OYsOckGXtOkrHnJJn5ek6s3GSxWCyWQKyRsFgsFksg5W4k7in0ARQh9pwkY89JMvacJDMvz0lZxyQsFovFEk65exIWi8ViCaEsjYSI3CAiO0Vkl4jcXujjKQQicoqI/FJEdojINhH5gLO9XUR+LiIvOI8LCn2sc42IVIrIZhH5kfN8lYg87pyTb4tITaGPca4RkTYR+Z6IPOdcM5eV+7UiIn/t3DtbReS/RaRuPl4rZWckRKQS+DfgRuAc4M0ick5hj6ogTAMfUkqdDVwK3Oqch9uBh5RSq4GHnOflxgeAHa7n/wx80jknJ4F3FuSoCsungZ8qpc4CLkSfn7K9VkSkC3g/sF4pdR5QCbyJeXitlJ2RAC4Bdiml9iilJoFvATcX+JjmHKXUEaXUU87vQ+ibvgt9Lu51drsX2FCYIywMIrIceCXwRee5AFcD33N2Kcdz0gK8FPgSgFJqUinVT5lfK+ilFupFpApoAI4wD6+VcjQSXcBB1/NDzrayRURWAuuAx4HFSqkjoA0JsKhwR1YQPgX8L2DWeb4Q6FdKTTvPy/F6OQ3oAb7iyHBfFJFGyvhaUUp1A/8CHEAbhwHgSebhtVKORkJ8tpVtipeINAH/A3xQKTVY6OMpJCLyKuC4UupJ92afXcvteqkCLgI+r5RaB4xQRtKSH0785WZgFbAMaERL2F5K/lopRyNxCDjF9Xw5cLhAx1JQRKQabSC+oZT6vrP5mIgsdf6+FDheqOMrAFcArxGRfWgZ8mq0Z9HmSApQntfLIeCQUupx5/n30EajnK+Va4C9SqkepdQU8H3gcubhtVKORuIJYLWThVCDDjbdX+BjmnMcrf1LwA6l1Cdcf7ofeLvz+9uBH8z1sRUKpdQdSqnlSqmV6OviYaXUW4BfAq93diurcwKglDoKHBSRNc6mVwDbKeNrBS0zXSoiDc69ZM7JvLtWyrKYTkRuQs8QK4EvK6XuLPAhzTkiciXwG+BZ4vr7R9Bxie8AK9A3whuUUicKcpAFRESuAv5WKfUqETkN7Vm0A5uBtyqlJgp5fHONiKxFB/NrgD3AO9CTzLK9VkTkH4A3ojMFNwN/gY5BzKtrpSyNhMVisViiUY5yk8VisVgiYo2ExWKxWAKxRsJisVgsgVgjYbFYLJZArJGwWCwWSyDWSFgsFoslEGskLJYQnJbqe0Wk3Xm+wHl+qs++K0XkT7L4rI9kc6wWSz6wRsJiCUEpdRD4PHCXs+ku4B6l1H6f3VcCGRsJdDGjxVJU2GI6iyUFTo+rJ4EvA+8C1jlt5r37PQacDexFt4n+DNqoXAXUAv+mlPpPp8/Rt4EWdPO896Dbk9+GroDf5rQDsVgKjjUSFksEROR64KfAdUqpnwfscxVOKw/n+S3AIqXU/xWRWuBR4A3A64A6pdSdziJYDUqpIREZVko1zcX/Y7FEpSr1LhaLBd0G+ghwHuBrJHy4DrhAREzDt1ZgNbrJ5JcdD2WjUmpLrg/WYskVNiZhsaTAaW53LXqZ17827bGjvBR4n1JqrfOzSin1M6XUr9ErvXUDXxeRt+XnyC2W7LFGwvL/2rtjlAZjMA7jz38Wr1G6uLg42jt4hg6OnsOpi2IdxUN4CVuL9AAdXTuJQhySQpGm0KHo8Py2fLyBfNNL3oQ32qO1gb6nPsq0Am6pL5LtsgZOt8YvwHXbMZBkkOSk3Yz6KKU8Utu1n7f4r02s9F+YJKT9xsBq6xziDhgmudwRuwC+k7wluaG21l4Cr0negQdqiXcEzJPMgCtg0uZPgUWS56P9jXQgD64lSV3uJCRJXd5ukg6U5Ax4+vX5s5Ry8RfrkY7JcpMkqctykySpyyQhSeoySUiSukwSkqQuk4QkqesHrleISyOaZzgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "f = X_test.dot(params['w']) + params['b']\n",
    "\n",
    "plt.scatter(range(X_test.shape[0]), y_test)\n",
    "plt.plot(f, color = 'darkorange')\n",
    "plt.xlabel('X_test')\n",
    "plt.ylabel('y_test')\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEKCAYAAADaa8itAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGnBJREFUeJzt3X+wXGWd5/H314REkF+BXJxMiAY0OuLUEvCKcVHLHXdDYEYDMzKCM5JCpuIqVunubC0wWkL5o0qdVUtKB0TNCi7DD0WKjKIxIqOlMwI3GPkhMokMQiSSQPg1otHAd/84zyXN5d4+py+nb98k71fVqT799HPOefp03/7c8zynT0dmIklSG54z6AZIknYfhookqTWGiiSpNYaKJKk1hookqTWGiiSpNYaKJKk1hookqTWGiiSpNTMH3YCpNnfu3Fy4cOGgmyFJu5R169Y9kJlDdfX2uFBZuHAhIyMjg26GJO1SIuIXTerZ/SVJao2hIklqjaEiSWqNoSJJao2hIklqjaEiSWqNoSJJao2h0tBnPgNXXDHoVkjS9GaoNHTBBfDVrw66FZI0vRkqkqTWGCo9yBx0CyRpejNUGooYdAskafozVCRJrTFUJEmtMVR64JiKJHVnqDTkmIok1TNUJEmtMVR6YPeXJHVnqDRk95ck1TNUJEmtMVQkSa0xVHrgmIokdWeoNOSYiiTVM1QkSa0xVHpg95ckdWeoNGT3lyTVM1QkSa0xVCRJrelbqETEgoi4PiLuiIjbI+I9pfy8iPhlRKwv0wkdy5wTERsj4s6IOK6jfFkp2xgRZ3eUHxYRN0TEhoi4IiJm9ev5gGMqklSnn0cqO4C/zcyXAUuAMyPiiPLYpzJzcZmuBSiPnQK8HFgG/ENEzIiIGcBngeOBI4BTO9bzsbKuRcBDwBn9ejKOqUhSvb6FSmZuzsyby/xjwB3A/C6LLAcuz8ztmfnvwEbgmDJtzMy7MvN3wOXA8ogI4E+Ar5blLwZO7M+zkSQ1MSVjKhGxEDgKuKEUvTsibomIVRExp5TNB+7tWGxTKZuo/GDg4czcMaa8b+z+kqTu+h4qEbEvcBXw3sx8FLgAeBGwGNgMfGK06jiL5yTKx2vDyogYiYiRrVu39vgMRtcxqcUkaY/S11CJiL2oAuXSzPwaQGben5lPZOaTwOepuregOtJY0LH4ocB9XcofAA6MiJljyp8hMy/KzOHMHB4aGmrnyUmSnqGfZ38F8EXgjsz8ZEf5vI5qJwG3lfnVwCkRMTsiDgMWATcCNwGLyples6gG81dnZgLXA28uy68ArunX85Ek1ZtZX2XSjgXeBtwaEetL2d9Rnb21mKqr6m7gHQCZeXtEXAn8lOrMsTMz8wmAiHg3sAaYAazKzNvL+s4CLo+IDwM/pgqxvnFMRZK661uoZOYPGH/c49ouy3wE+Mg45deOt1xm3sXO7rO+ckxFkur5jXpJUmsMlR7Y/SVJ3RkqDdn9JUn1DBVJUmsMFUlSawyVHjimIkndGSoNOaYiSfUMFUlSawyVHtj9JUndGSoN2f0lSfUMFUlSawwVSVJrDJUeOKYiSd0ZKg05piJJ9QwVSVJrDJUe2P0lSd0ZKg3Z/SVJ9QwVSVJrDBVJUmsMlR44piJJ3RkqDTmmIkn1DBVJUmsMlR7Y/SVJ3RkqDdn9JUn1DBVJUmsMFUlSawyVHjimIkndGSoNOaYiSfUMFUlSawyVHtj9JUndGSoN2f0lSfX6FioRsSAiro+IOyLi9oh4Tyk/KCLWRsSGcjunlEdEnB8RGyPilog4umNdK0r9DRGxoqP8FRFxa1nm/Ag/+iVpkPp5pLID+NvMfBmwBDgzIo4Azgauy8xFwHXlPsDxwKIyrQQugCqEgHOBVwHHAOeOBlGps7JjuWV9fD6SpBp9C5XM3JyZN5f5x4A7gPnAcuDiUu1i4MQyvxy4JCs/Ag6MiHnAccDazNyWmQ8Ba4Fl5bH9M/NfMzOBSzrW1afn1M+1S9Kub0rGVCJiIXAUcAPw/MzcDFXwAIeUavOBezsW21TKupVvGqd8vO2vjIiRiBjZunXrJJ/DpBaTpD1K30MlIvYFrgLem5mPdqs6TllOovyZhZkXZeZwZg4PDQ3VNVmSNEl9DZWI2IsqUC7NzK+V4vtL1xXldksp3wQs6Fj8UOC+mvJDxynvG7u/JKm7fp79FcAXgTsy85MdD60GRs/gWgFc01F+WjkLbAnwSOkeWwMsjYg5ZYB+KbCmPPZYRCwp2zqtY119eD79WrMk7T5m9nHdxwJvA26NiPWl7O+AjwJXRsQZwD3AyeWxa4ETgI3A48DpAJm5LSI+BNxU6n0wM7eV+XcCXwL2Br5ZJknSgPQtVDLzB4w/7gHwhnHqJ3DmBOtaBawap3wE+ONn0UxJUov8Rn0PHFORpO4MlYYcU5GkeoaKJKk1hkoP7P6SpO4MlYbs/pKkeoaKJKk1hookqTWGSg8cU5Gk7gyVhhxTkaR6hookqTWGSg/s/pKk7gyVhuz+kqR6hookqTWGiiSpNYZKDxxTkaTuDJWGHFORpHqGiiSpNYZKD+z+kqTuDJWG7P6SpHqGiiSpNYaKJKk1hkoPHFORpO4MlYYcU5Gkeo1CJSLeExH7R+WLEXFzRCztd+MkSbuWpkcqb8/MR4GlwBBwOvDRvrVqmrL7S5K6axoqo50/JwD/NzN/0lG2R7D7S5LqNQ2VdRHxbapQWRMR+wFP9q9ZkqRd0cyG9c4AFgN3ZebjEXEQVReYJElPaXqk8mrgzsx8OCL+Gng/8Ej/mjU9OaYiSd01DZULgMcj4kjgfwO/AC7pW6umIcdUJKle01DZkZkJLAc+nZmfBvbrtkBErIqILRFxW0fZeRHxy4hYX6YTOh47JyI2RsSdEXFcR/myUrYxIs7uKD8sIm6IiA0RcUVEzGr6pCVJ/dE0VB6LiHOAtwHfiIgZwF41y3wJWDZO+acyc3GZrgWIiCOAU4CXl2X+ISJmlO18FjgeOAI4tdQF+FhZ1yLgIapxn76y+0uSumsaKm8BtlN9X+VXwHzg77stkJnfB7Y1XP9y4PLM3J6Z/w5sBI4p08bMvCszfwdcDiyPiAD+BPhqWf5i4MSG25oUu78kqV6jUClBcilwQET8GfDbzJzsmMq7I+KW0j02p5TNB+7tqLOplE1UfjDwcGbuGFM+rohYGREjETGydevWSTZbklSn6WVa/hK4ETgZ+Evghoh48yS2dwHwIqrTkzcDnxjdxDh1cxLl48rMizJzODOHh4aGemuxJKmxpt9TeR/wyszcAhARQ8B32Nn91Ehm3j86HxGfB75e7m4CFnRUPRS4r8yPV/4AcGBEzCxHK531+8YxFUnqrumYynNGA6V4sIdlnxIR8zrungSMnhm2GjglImZHxGHAIqojo5uAReVMr1lUg/mry5lo1wOjR0srgGt6bU9vbe/n2iVp99D0SOVbEbEGuKzcfwtwbbcFIuIy4PXA3IjYBJwLvD4iFlN1Vd0NvAMgM2+PiCuBnwI7gDMz84mynncDa4AZwKrMvL1s4izg8oj4MPBj4IsNn4skqU8iG/bpRMRfAMdSjWd8PzOv7mfD+mV4eDhHRkZ6Xu744+HBB+HGG/vQKEma5iJiXWYO19VreqRCZl4FXPWsWrULs/tLkup1DZWIeIzxz6oKIDNz/760SpK0S+oaKpnZ9VIskiR18jfqe+ApxZLUnaHSkGMqklTPUJEktcZQ6YHdX5LUnaHSkN1fklTPUJEktcZQkSS1xlDpgWMqktSdodKQYyqSVM9QkSS1xlDpgd1fktSdodKQ3V+SVM9QkSS1xlCRJLXGUOmBYyqS1J2h0pBjKpJUz1CRJLXGUOmB3V+S1J2h0pDdX5JUz1CRJLXGUJEktcZQ6YFjKpLUnaHSkGMqklTPUJEktcZQ6YHdX5LUnaHSkN1fklTPUJEktaZvoRIRqyJiS0Tc1lF2UESsjYgN5XZOKY+IOD8iNkbELRFxdMcyK0r9DRGxoqP8FRFxa1nm/AiPJSRp0Pp5pPIlYNmYsrOB6zJzEXBduQ9wPLCoTCuBC6AKIeBc4FXAMcC5o0FU6qzsWG7stlrnmIokdde3UMnM7wPbxhQvBy4u8xcDJ3aUX5KVHwEHRsQ84DhgbWZuy8yHgLXAsvLY/pn5r5mZwCUd6+oLj4Mkqd5Uj6k8PzM3A5TbQ0r5fODejnqbSlm38k3jlI8rIlZGxEhEjGzduvVZPwlJ0vimy0D9eMcBOYnycWXmRZk5nJnDQ0NDk2yi3V+SVGeqQ+X+0nVFud1SyjcBCzrqHQrcV1N+6DjlfWP3lyTVm+pQWQ2MnsG1Arimo/y0chbYEuCR0j22BlgaEXPKAP1SYE157LGIWFLO+jqtY12SpAGZ2a8VR8RlwOuBuRGxieosro8CV0bEGcA9wMml+rXACcBG4HHgdIDM3BYRHwJuKvU+mJmjg//vpDrDbG/gm2WSJA1Q30IlM0+d4KE3jFM3gTMnWM8qYNU45SPAHz+bNvbKMRVJ6m66DNRPe46pSFI9Q0WS1BpDpQd2f0lSd4ZKQ3Z/SVI9Q0WS1BpDRZLUGkOlB46pSFJ3hkpDjqlIUj1DRZLUGkOlB3Z/SVJ3hkpDdn9JUj1DRZLUGkNFktQaQ6WhCHjyyUG3QpKmN0OloRkz4IknBt0KSZreDJWGDBVJqmeoNGSoSFI9Q6UhQ0WS6hkqDRkqklTPUGnIUJGkeoZKQzNmwI4dg26FJE1vhkpDHqlIUj1DpSFDRZLqGSoNGSqSVM9QaWjmTENFkuoYKg15pCJJ9QyVhmbMqG69qKQkTcxQaWg0VDxakaSJGSoN7b13dfv444NthyRNZ4ZKQ/PnV7f33jvYdkjSdDaQUImIuyPi1ohYHxEjpeygiFgbERvK7ZxSHhFxfkRsjIhbIuLojvWsKPU3RMSKfrZ5eLi6/fKX+7kVSdq1DfJI5b9k5uLMLB/XnA1cl5mLgOvKfYDjgUVlWglcAFUIAecCrwKOAc4dDaJ+ePGLYcUK+NSn4I47+rUVSdq1Tafur+XAxWX+YuDEjvJLsvIj4MCImAccB6zNzG2Z+RCwFljWzwZ+/OPwvOfBu97lWWCSNJ5BhUoC346IdRGxspQ9PzM3A5TbQ0r5fKBzJGNTKZuo/BkiYmVEjETEyNatWyfd6EMOqYLln/8ZPv3pSa9GknZbgwqVYzPzaKqurTMj4nVd6sY4Zdml/JmFmRdl5nBmDg8NDfXe2g5/8zewfDmcdRbcfPOzWpUk7XYGEiqZeV+53QJcTTUmcn/p1qLcbinVNwELOhY/FLivS3lfRcAXvgBDQ3DSSfCrX/V7i5K065jyUImI50XEfqPzwFLgNmA1MHoG1wrgmjK/GjitnAW2BHikdI+tAZZGxJwyQL+0lPXd3LmwejU88AC88Y3w619PxVYlafqbOYBtPh+4OiJGt/+PmfmtiLgJuDIizgDuAU4u9a8FTgA2Ao8DpwNk5raI+BBwU6n3wczcNlVP4hWvgMsuq45Wli2Db3wD9t9/qrYuSdNTZI47DLHbGh4ezpGRkdbW95WvwFvfCkcdVR29/MEftLZqSZo2ImJdx1dAJjSdTineJZ18Mlx1Fdx+e3X0csMNg26RJA2OodKCN70J/uVfYNYseM1r4AMfgN/9btCtkqSpZ6i05MgjYd06OPVU+NCHYPFi+Kd/gj2sd1HSHs5QadFBB8Ell8DXvw47dlRHMK99LVxzjZfMl7RnMFT64E//tBpjufBCuOceOPFEeMlL4CMfgZ//fNCtk6T+MVT6ZK+94B3vgLvuqs4QO/RQeP/7qwtTvvKVcN558MMfwu9/P+iWSlJ7PKV4Ct1zD1x5ZRUyN91Ujbfstx8sWVIFzej0h39YfXNfkqaLpqcUGyoDsm0bfPe78J3vVKch33rrznGX/feHP/qjanrpS2HRInjBC2DBgup7MM/x+FLSFDNUJjBdQmWs3/wG1q+vziD72c92Tr/85dPrzZxZdaUtWADz5lXXIBudDjlk5/zcuXDAATB79mCej6TdS9NQGcRlWjSOvfeGV7+6mjo99lg1LnPvvTune+6pbtevh61b4aGHJl7v7NlVuEw07bsv7LNP9Tsx++yzc+p2f9Ysj5Ykjc9Qmeb226/6DsyRR05c5/e/hwcfhC1bqpDZurW6/8gjO6eHH945v3nzzvlf/3py36WZObMKrNmzq5DpNj+2bNasavmZM2HGjJ3zE01N68yYUYVd3RTRrF4v0+gYWERv89LuxlDZDey1VzXWMpnrjmXC9u1VuDz++M5p7P3Osu3bqysGbN8+8fzo7UMPPf3x7dur7/CMnZ54orrdE9WFT9OQmkywTTb0pvPj07ltg378xz/uf5e4obKHi4DnPreaDj54sG3JrH6muVvodJuefLL5NLqtNqbREywydx711c0Pql4vy3R7nabr49O5bdPh8ak4QjZUNG1E7OzG8gQDadfkcKskqTWGiiSpNYaKJKk1hookqTWGiiSpNYaKJKk1hookqTWGiiSpNXvcVYojYivwi0kuPhd4oMXmtMV29cZ29cZ29WZ3bdcLM3OortIeFyrPRkSMNLn081SzXb2xXb2xXb3Z09tl95ckqTWGiiSpNYZKby4adAMmYLt6Y7t6Y7t6s0e3yzEVSVJrPFKRJLUnM51qJmAZcCewETi7T9tYAFwP3AHcDrynlJ8H/BJYX6YTOpY5p7TpTuC4uvYChwE3ABuAK4BZDdt2N3Br2f5IKTsIWFvWtRaYU8oDOL9s+xbg6I71rCj1NwArOspfUda/sSwbDdr00o59sh54FHjvIPYXsArYAtzWUdb3/TPRNmra9ffAz8q2rwYOLOULgd907LcLJ7v9bs+xS7v6/roBs8v9jeXxhQ3adUVHm+4G1g9gf0302TDw99i4fw/9+IDcnSZgBvBz4HBgFvAT4Ig+bGfe6IsP7Af8G3BE+WP7X+PUP6K0ZXb5I/p5aeuE7QWuBE4p8xcC72zYtruBuWPKPj76hwycDXyszJ8AfLO8sZcAN3S8Oe8qt3PK/OgfwY3Aq8sy3wSOn8Rr9CvghYPYX8DrgKN5+odR3/fPRNuoaddSYGaZ/1hHuxZ21huznp62P9FzrGlX31834F2UD3/gFOCKunaNefwTwAcGsL8m+mwY+Hts3Off64ffnjaVHb2m4/45wDlTsN1rgP/W5Y/tae0A1pS2jtve8mZ5gJ0fKE+rV9OWu3lmqNwJzCvz84A7y/zngFPH1gNOBT7XUf65UjYP+FlH+dPqNWzfUuCHZX4g+4sxHzJTsX8m2ka3do157CTg0m71JrP9iZ5jzf7q++s2umyZn1nqRbd2dZQHcC+waBD7a8w2Rj8bpsV7bOzkmEq9+VRvplGbSlnfRMRC4CiqQ3SAd0fELRGxKiLm1LRrovKDgYczc8eY8iYS+HZErIuIlaXs+Zm5GaDcHjLJds0v82PLe3EKcFnH/UHvL5ia/TPRNpp6O9V/paMOi4gfR8T3IuK1He3tdfuT/Zvp9+v21DLl8UdK/SZeC9yfmRs6yqZ8f435bJiW7zFDpV6MU5Z921jEvsBVwHsz81HgAuBFwGJgM9UheLd29VrexLGZeTRwPHBmRLyuS92pbBcRMQt4E/CVUjQd9lc306IdEfE+YAdwaSnaDLwgM48C/ifwjxGx/yS3P5llpuJ1ezb78lSe/o/LlO+vcT4bel3flLzHDJV6m6gGykYdCtzXjw1FxF5Ub5pLM/NrAJl5f2Y+kZlPAp8Hjqlp10TlDwAHRsTMXp9HZt5XbrdQDe4eA9wfEfNKu+dRDXBOpl2byvzY8qaOB27OzPtLGwe+v4qp2D8TbaOriFgB/BnwV1n6NTJze2Y+WObXUY1XvGSS2+/5b2aKXrenlimPHwBs69aujrp/TjVoP9reKd1f4302TGJ9U/IeM1Tq3QQsiojDyn/FpwCr295IRATwReCOzPxkR/m8jmonAbeV+dXAKRExOyIOAxZRDbaN297y4XE98Oay/Aqqvtm6dj0vIvYbnacav7itbH/FOOtaDZwWlSXAI+WweQ2wNCLmlK6NpVR93ZuBxyJiSdkHpzVpV4en/Qc56P3VYSr2z0TbmFBELAPOAt6UmY93lA9FxIwyf3jZP3dNcvsTPcdu7ZqK162zvW8GvjsaqjX+K9WYw1NdRFO5vyb6bJjE+qbkPdbXwebdZaI6m+LfqP4beV+ftvEaqkPOW+g4rRL4MtWpfreUF3hexzLvK226k44zpiZqL9WZMjdSnTb4FWB2g3YdTnVmzU+oTmd8Xyk/GLiO6lTD64CDSnkAny3bvhUY7ljX28u2NwKnd5QPU32I/Bz4DA1OKS7L7QM8CBzQUTbl+4sq1DYDv6f6r++Mqdg/E22jpl0bqfrVn3YqLPAX5fX9CXAz8MbJbr/bc+zSrr6/bsBzy/2N5fHD69pVyr8E/Pcxdadyf0302TDw99h4k9+olyS1xu4vSVJrDBVJUmsMFUlSawwVSVJrDBVJUmsMFWmai4jXR8TXB90OqQlDRZLUGkNFaklE/HVE3BgR6yPicxExIyL+IyI+ERE3R8R1ETFU6i6OiB9FdQHFq8s3nImIF0fEdyLiJ2WZF5XV7xsRX42In0XEpeWbz0TERyPip2U9/2dAT116iqEitSAiXga8herim4uBJ4C/Ap5HdW2yo4HvAeeWRS4BzsrM/0T1refR8kuBz2bmkcB/pvqGN1RXpn0v1e9oHA4cGxEHUV3S5OVlPR/u77OU6hkqUjveQPXreTdFxPpy/3DgSXZeiPD/Aa+JiAOofnHxe6X8YuB15Rpr8zPzaoDM/G3uvD7XjZm5KasLLq6n+j2PR4HfAl+IiD8HnrqWlzQohorUjgAuzszFZXppZp43Tr1u10Ua7xLko7Z3zD9B9SNUO6iu5nsVcCLwrR7bLLXOUJHacR3w5og4BCAiDoqIF1L9jY1eMfetwA8y8xHgodj5w05vA76X1W9kbIqIE8s6ZkfEPhNtMKrf1zggM6+l6hpb3I8nJvViZn0VSXUy86cR8X6qX8h8DtWVbs8Efg28PCLWUf3S4FvKIiuAC0to3AWcXsrfBnwuIj5Y1nFyl83uB1wTEc+lOsr5Hy0/LalnXqVY6qOI+I/M3HfQ7ZCmit1fkqTWeKQiSWqNRyqSpNYYKpKk1hgqkqTWGCqSpNYYKpKk1hgqkqTW/H+XISJuumqCdQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(loss_his, color = 'blue')\n",
    "plt.xlabel('epochs')\n",
    "plt.ylabel('loss')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(442, 11)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.utils import shuffle\n",
    "X, y = shuffle(data, target, random_state=13)\n",
    "X = X.astype(np.float32)\n",
    "data = np.concatenate((X, y.reshape((-1,1))), axis=1)\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(353, 10) (353, 1) (89, 10) (89, 1)\n",
      "epoch 10000 loss 5667.323261\n",
      "epoch 20000 loss 5365.061628\n",
      "epoch 30000 loss 5104.951009\n",
      "epoch 40000 loss 4880.559191\n",
      "epoch 50000 loss 4686.465012\n",
      "epoch 60000 loss 4518.097942\n",
      "epoch 70000 loss 4371.603217\n",
      "epoch 80000 loss 4243.728424\n",
      "epoch 90000 loss 4131.728130\n",
      "five kold cross validation score is 4033.2928969313452\n",
      "valid score is 3558.147428823725\n",
      "(353, 10) (353, 1) (89, 10) (89, 1)\n",
      "epoch 10000 loss 5517.461607\n",
      "epoch 20000 loss 5181.523292\n",
      "epoch 30000 loss 4897.000442\n",
      "epoch 40000 loss 4655.280406\n",
      "epoch 50000 loss 4449.235361\n",
      "epoch 60000 loss 4272.963989\n",
      "epoch 70000 loss 4121.578253\n",
      "epoch 80000 loss 3991.027359\n",
      "epoch 90000 loss 3877.952427\n",
      "five kold cross validation score is 3779.575673476155\n",
      "valid score is 4389.749894361285\n",
      "(354, 10) (354, 1) (88, 10) (88, 1)\n",
      "epoch 10000 loss 5677.526932\n",
      "epoch 20000 loss 5331.464606\n",
      "epoch 30000 loss 5038.541403\n",
      "epoch 40000 loss 4789.859079\n",
      "epoch 50000 loss 4578.047735\n",
      "epoch 60000 loss 4397.001364\n",
      "epoch 70000 loss 4241.659289\n",
      "epoch 80000 loss 4107.825490\n",
      "epoch 90000 loss 3992.019241\n",
      "five kold cross validation score is 3891.3610347157187\n",
      "valid score is 3844.049401675179\n",
      "(354, 10) (354, 1) (88, 10) (88, 1)\n",
      "epoch 10000 loss 5366.435738\n",
      "epoch 20000 loss 5117.050627\n",
      "epoch 30000 loss 4902.593420\n",
      "epoch 40000 loss 4717.676380\n",
      "epoch 50000 loss 4557.767422\n",
      "epoch 60000 loss 4419.053021\n",
      "epoch 70000 loss 4298.323163\n",
      "epoch 80000 loss 4192.874781\n",
      "epoch 90000 loss 4100.430691\n",
      "five kold cross validation score is 4019.0791769506664\n",
      "valid score is 4300.006957377207\n",
      "(354, 10) (354, 1) (88, 10) (88, 1)\n",
      "epoch 10000 loss 5592.024970\n",
      "epoch 20000 loss 5278.867272\n",
      "epoch 30000 loss 5010.647792\n",
      "epoch 40000 loss 4780.303445\n",
      "epoch 50000 loss 4581.916469\n",
      "epoch 60000 loss 4410.526730\n",
      "epoch 70000 loss 4261.974936\n",
      "epoch 80000 loss 4132.771612\n",
      "epoch 90000 loss 4019.987616\n",
      "five kold cross validation score is 3921.1718979156462\n",
      "valid score is 3909.701727644273\n"
     ]
    }
   ],
   "source": [
    "from random import shuffle\n",
    "\n",
    "def k_fold_cross_validation(items, k, randomize=True):\n",
    "    if randomize:\n",
    "        items = list(items)\n",
    "        shuffle(items)\n",
    "\n",
    "    slices = [items[i::k] for i in range(k)]\n",
    "\n",
    "    for i in range(k):\n",
    "        validation = slices[i]\n",
    "        training = [item\n",
    "                    for s in slices if s is not validation\n",
    "                    for item in s]\n",
    "        training = np.array(training)\n",
    "        validation = np.array(validation)\n",
    "        yield training, validation\n",
    "\n",
    "\n",
    "for training, validation in k_fold_cross_validation(data, 5): \n",
    "    X_train = training[:, :10]\n",
    "    y_train = training[:, -1].reshape((-1,1))\n",
    "    X_valid = validation[:, :10]\n",
    "    y_valid = validation[:, -1].reshape((-1,1))\n",
    "    loss5 = []\n",
    "    #print(X_train.shape, y_train.shape, X_valid.shape, y_valid.shape)\n",
    "    loss, params, grads = linar_train(X_train, y_train, 0.001, 100000)\n",
    "    loss5.append(loss)\n",
    "    score = np.mean(loss5)\n",
    "    print('five kold cross validation score is', score)\n",
    "    y_pred = predict(X_valid, params)\n",
    "    valid_score = np.sum(((y_pred-y_valid)**2))/len(X_valid)\n",
    "    print('valid score is', valid_score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(353, 10) (353, 1) (89, 10) (89, 1)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.datasets import load_diabetes\n",
    "from sklearn.utils import shuffle\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "diabetes = load_diabetes()\n",
    "data = diabetes.data\n",
    "target = diabetes.target \n",
    "X, y = shuffle(data, target, random_state=13)\n",
    "X = X.astype(np.float32)\n",
    "y = y.reshape((-1, 1))\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)\n",
    "print(X_train.shape, y_train.shape, X_test.shape, y_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficients: \n",
      " [[   7.00153867 -244.27466725  488.58606379  301.83242902 -693.67792807\n",
      "   364.02013146  106.15852423  289.24926974  645.158344     50.77526251]]\n",
      "Mean squared error: 3371.88\n",
      "Variance score: 0.54\n",
      "0.5392080506325068\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADuCAYAAAAOR30qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXmYXFWZ/7/Vne5OpzvpTtKdhBC6m4RAguwEWQKyiUSFQeFBGMI4iLgAIzrzU+c3k/nhCqMwow/KICKjoLQ6jgoCgqyRPQKRnQQNWZqsZO/O1uv9/XH61nnPueeudbeqej/Pk6dvVVdu3a6693u/533f856CZVlgGIZhsqcm6wNgGIZhBCzIDMMwOYEFmWEYJiewIDMMw+QEFmSGYZicwILMMAyTE1iQGYZhcgILMsMwTE5gQWYYhskJY8K8uK2tzerq6kroUBiGYSqTpUuXbrEsq93vdaEEuaurCy+++GL0o2IYhqlCCoXCmiCv45AFwzBMTmBBZhiGyQksyAzDMDmBBZlhGCYnsCAzDMPkBBbkcqK7G+jqAmpqxM/u7qyPiGGYGAlV9sZkSHc38OlPA3v2iMdr1ojHALBwYXbHxTBMbLBDLhcWLZJibLNnj3ieYZiKgAW5XOjpCfc8wzBlBwtyudDREe55hmHKDhbkcuG664Bx49Tnxo0TzzMMUxGwIJcLCxcCt90GdHYChYL4edtt6SX0uMKDYRKHqyzKiYULs6mo4AoPhkkFdsiMP1zhER4eUTARYIfM+MMVHuHgEQUTEXbIjD9c4REOHlEwEWFBZvzhCo9w8IiCiQgLMuNP1hUe5QaPKJiIsCAzwVi4EFi9GhgZET9ZjN3hEQUTERZkhokbHlEwEeEqC4ZJgqxqxpmyhh0ywzBMTmBBZhiGyQksyADPqmIYJhdwDJlnVTEMkxPYIfOsKoZhcgILMs+qYhgmJ7Ag86wqhmFyAgsyz6rKDk6mMowCCzLPqsoGO5m6Zg1gWTKZyqLMVDEFy7ICv3jevHnWiy++mODhMFVDV5cQYZ3OTtErg2EqiEKhsNSyrHl+r2OHzGQDJ1MZxgELMpMNnExlGAcsyEw2cDKVYRywIDPZwMlUhnHAgsxkR9Sm91wux1QoLMh5g8XGGy6XYyoYFuQ8wWLjD/ceYSoYFuQ8kZTYVJLr5nI5poJhQc4TSYhNpbluLpdjKhgW5DyRhNhU2hCfy+WYCoYFOU8kITaVNsTncjmmgmFBzhNJiE0lDvGjlssxTM5hQc4bcYsND/EZpmxgQa50eIjPMGUDC3I1wEP8XLNkCXDllcCzz6bwZpVUAlmB8KrTDJMhlgVccgmwahXwwAPiflkoJPRmvMJ67mGHzDAZ0t8vxBgQhS/9/Qm+WaWVQFYgLMgMkyF9ferjRAU5bAkkhzdShwU5TvgEZkKiC/K+fQm+WZgSyEqb4VkmsCDHBZ/ATARSdchhSiArIbxRhgaJBTkuKuEEZlKnt1d9nKgghymBLPcZnmVqkFiQ46LcT2AmExwO+dQPJCsaQUsgy32GZ5kaJBbkuCj3E5jJhL4Hn1Ye92/Ymg8nV+4zPMvUILEgx0W5n8BMJvT9/D7lcT8a8uHkyn2GZ5kaJBbkuCj3E5jJhL5tg8rjfRgrNvycXBoJq3Ke4VmmBim/glyGGdKyPoGZTOhrmaE87keD2PBycmWasEqVMjVI+RRkPuFcefttYP36rI+CiYu++QuUx/1ocHdytkm59NKyTFilThkapHwKcplmSJPmkUeAgw4S5mnZsqyPhomD3hmHKo/722aYnRw1KW7kPGHF+JPP5kJlmiFNmnvuET+Hh0Ujmrlzsz0epnQcZW//eTNgMnImk6KT84QV408+HXKZZkiTZudOub1rV3bHwcRH4KnTfmakDBJWjD/5FOQyzZAmDRVkP7NULlhW1keQLYGnTnuZkTJJWDH+5FOQyzRDmjR0mu3u3dkdR1w8+SSw//7AGWcAg4P+r69EAguym0m5666ySVgx/uRTkIGyzJAmDXXInoJcJiWDN90EbNgALF4MPPxw1keTDYEFmU1KVZDPpB5jJJBDLqNVITZsMG9XE6G6vS1cmLvvkImX/DpkxkGgGHIZlQxu2ya3t2zJ7jiyRO/2lmg/ZCb3VJUg//KXwHvfC/z4x1kfSXgsK6BDLqOSwa1b5XY1CvLwsPPemWj7TSb3VJUgX3MN8MILwFVXOYeKeWfvXmBoSD52FeQyKRkcGVEdMhXn1Mko5m4qXWRBrm6qRpD7+4HNm+X2Cy9kezxhoeEKwEOQEyoZtCzg3nuBn/9cvTFEpbdXiLJNZg45w2n6JlOQtCCPjABXXAGcdhrP9swjVSPI+sn/3AU35rb6wIQea3SNISeUjV+8GDjvPLGbn/60pF0BUN0xkKEgZxhzz0KQH38c+O//Bp54AvjWt5J9LyY8VSPIvT/7nfJ4yY45ZdWwKLBDBhIpGXzmGbn95JMl784RoshMkDOMuZsEOemk3rp1cnvFCu/Xbt4M3HknN7NKk+oR5Bt/qDxeghNg5bT6wITukNOeGGKHewD1oo6K7pAziyFnGHPPwiHTuLWf0H7sY8BllwHvfz/PqEyL6hHkDaqCbUE7VuAgETPM8eQJG90h79mjxmCTJm5B1gV4+3aP2HSSSbcMp+nrN1kgXUHesMFdaIeGRFgDELHmHTuSPa5ckeHEquoR5PZZjueW4ASxUQb9lk0Xb5o1q1SQ164tfX+6Q3Z7LvGkW4Yz4LJ2yP397kK7ZYsq1nv3JntcuSHjXuzVI8gXftLx3HM4UT7IefhCd8hAumELGuPt6yu9bNAUojCGLdJIumU0TT9rQQaA9bfdb3SDmzapr6saQc54YlX1CPLh8x3PKYIMwFpTWiJn1SrgyivjqULQyVqQqUMGSg9bmNywMbFXRhNdwpIHQd7wlVuNbrBqBTnj8616BNkw5H8VR2A3xsECcCVuQVNhN268Mfp7LFoE3HqrSISsWhV9PyZMx5+WIFtW/IJscsNGQS6TiS5RyKLKwuGQ+yepT4y6waoV5IzPt6oW5BHU4gUch9/gAtyKK7HXasQ3vhE9o7x8ufhpWcBf/xr9WE2YHHJaPZF7e53tMVNzyBXcGzsXDhn7OV/U01O9gpzx+VY1gkwFbQzpcfcozsI/1d5UfNzXF70Ei4q+6QZQClk6ZJNQJuGQjZ97BbedzEKQ9XNmPaY7X9TRUb2CnPH5VjXtN6mgnXgi8NRTYvuGukUO97dmDdDWFv496AUWtyBnGUPWwxVAig4ZqNi2k2EF2bLElP+2NmDmzGjv6XDItTOAYfLEqBvc+Af1dZWyQk0gMjzfqsYhU4E8+2y5bVqpImr8nl5gJgEthSwdchKCHDiGXMGErUPu7gaOPx6YPTt6SMwRQ551itENVq1DzpiqdMjHHQe0trrXYHqttO7G0JB60qbhkNNyLXEL8vCw+bOvNkEOm9R76CHxc2REbM+eHf49HQ55aIoo9dNgQc6GqnTIra3CadgUCmKaqE0UQdZP9LRDFn19wH/8B/CrX8X7voBZkEuZHLJjhzlxmmkLzgwwCfLIiPuMRXqOvftutPd0CLLLbL2KE+QyWdasKgV5wgTg9NPl489+FvjoR+XjKCELXYDTTup997vAl74EXHQR8OKL8b63yblu2hS9DSeNHxcK3u9TybhNrnELWyQhyHv3Om/2w8POm3BZC3LGs+/CULWCfNVVwGc+A1x9NXDjjSJ8ZhPFIesXV5wx5OFhczNzGrJ46SW5vXRpfO8NmB3yyAiwcWO0/VEnTD93N0EeHBQ3zPe8B/jzn6O9Zx5JW5AHBsw5E309w61bnX1SylqQy2hZs6oV5PHjxSSOm28GmprUuu84BDlOh+x24VKHTGOycQ/9TYIMRI8jU4c8a5Z0yTt2mF334sXAPfcAb74pRgGVgGWp3+vEiXI7KUE23dQBZ9c3PVwBlLkgl9Fsz/wIcoIxnqEheYMsFIQA6+y3H1BXJ7a3bAmfMEtSkN3cNhVk+pq4h/5UkCdPlttRBZneMNrbgUlkspipHI6+/xNPiM5w5cDgoEvDJAjRtW8+dXXCJNi4JfaSEmTdIVecIJcy+y7l2HM+BDnhGA8VywkT1LilTU0NcMAB8nHYm2dcMeShIeCWW4CbbpLDS7d9uQlykg75qKPkdhwOefJktebbdDOhYjA8DDz4YLT3TZMtW8T5NH26rI6g0O90/HigoUE+DuKQTaLpRykOuazrkKPOvssg9pwPQU44xqOHK9woJWwRl0O+5x4R1/7CF2STIjeHTD+yJEMWVCTjEGR6fJMmqa7bT5AB4He/c74mb/z610LY+vuBn/3M+Xt6vkQR5N7e8H0vqtYhR519l0HsOR+CnHCMJ6gg0wRT2LeOK6n36qty+09/Ej/p8VN3bztky0rOIe/bJy/kujpg7lz5uyQcsunYdTF48EGRoMozf/mL3DbVXOujNj9BHhpyCrBbbN+NoIJsStZmKshxhA2itFjNIPacD0FOuMNSFEEu1SH39UVb0YNeNO+8I35SsW1vl9u2IO/ZI4byNnEKMr3o29qAGTPk47gccpiQBSA+2z/+Mdp7pwWdSWcaLYV1yKZZmWHjyPTcoiP4XCf1sixZy6DzWz4EOeEOS2mELPSLzrKiTW2mF409+YIK8nTSC8bev+7GIwuywYlQQW5vB/bf33l8YdEdsl/IwhS/vPfeaO+dFtQhBxHksWPlY1MowuRuSxFkOssv1yGLLEvWMuj8lg9BTrjDUhYhC/19o+zHdsh0P/uRbon2eaoPibdvVx1zIFycyJZfPVZ8iS7I69ZFa1VKbxhRQhaAEOS8Lrw5NASsXCkfx+GQTYIcNrHnJsjr16ufZa4EOcuStQw6v+VDkIFkltEZdXy9F3+6+FRaIQsgWhyZXjQ7d4r90v1QQXZzyJYVYVFKFyey+XZpRdvaxLRz2zTs2RPtb6QOOUrIAhA3q5dfDv/eabBmjVpPnZQgl+KQ998faGwU23v2qMeTK0HOeoGClJf3yo8gxw1xfL0YX3x6wrplrv+Flr2tXRtuanBcDlm/8NauVfcTJGQBRAhbuDiOzVvlKdLeLoyC7pLDojvkMFUWtogA+a220Dux9fY63XzYsre4BXn8ePXmbseRR0bM+81MkBMKG6xZA3zwg8Dll5tnL2ZF5QoycXy9kLZ4wnOGotBRxo4Fpk4V28PDztiaFybxjUOQ33nH3yGb3HBoQXZxHJsnyNW67YRiKYI8OCg/l0IBaGkJ55DPOUdulxRHDpq5j5Dhp/FjQPzNusiGrbKIQ5BpTqO5Wb252+f6tm3mcFfgOuS4J1IkFDa44QbgD38AfvIT4L77SjvEOKlcQSaOjwpyy07v2FPUsEVaDtm+YQAi+TMyEpNDdnEim4+VzaNLFuTubuyYeUzx4cSmftTWhoshf/SjckblSy9FmyAROHMfMcNv6lWsnwtZJ/Wam80OmX6edDp3IIecVEVEAmEDWl66YkXJu4uNyhVk4vgUhzzJuwV01EqLJGLIgNMhT5yoDtvd4rihBdnFiWxpldmfkgR59GLdulZarUm73wG6u0OFLKZMAebMkY/1kq1ABM3cR8zwRxHktEMWTU1mh0wFuatLbgcS5DJq4vP223I7apOsJKhcQSaObydaik9PWHiu53+LWmmRRJUFIBwyFdwJE9ReHLt3xxSyAIxORK9DBiIK8ujFug2yccVkawuwaBEmTvRuMKTHkKlzi9TXImjmPmKGPy1BLqXKQnfItiBTcdIF2beqpUya+OzerYYjWZDTgDg+xSF/+BTP/5ZlyMKyzA6Z7qelxSnIfg55wwbgi1+MNnLU65CBiLXIoxflVkg7PAnbgJ4e1NZ6NxiKXZCDZu4jZPgHBowLcCTmkMOU/oUNWUyfLsNDlhVgdmTWFREBoSWJAAtyeow6vt7jP1B8yqvsDYgWstBbKdqEFeT+fmdCxeSQaag3SMhi0SLgP/8TuPRStW9yEEyCHGm23ugHqzhkbC0+75XY8xJkt25qngTN3EfI8K9caZ6hmYQgDw2FK2/UBdkvZDFtmhoe8w1bZDCRIgo0XAGwIKdO0IkhQLSQxe7dZqcSNoZsuuiCOGS/kAVdQeThb78UOAs+NCQFr1CQ5WmRQhajF6vikMf0FS9WrzgyFYJx42JwyEEz9xEy/G6Lj3p1AwxSZeHWEztMHDmsQ546NaQgZzCRIgp6Ei9PglwVi5xGFWQ7WWxq10lxu1jCOmSTINN91NeLbHzYkIU94w8Anvn1emB41PrbWXDAeNHoEzhqa8X2tGniM7EsIQhDQ8AYvzNpdP/brtoGjP5Nk8+bDyw8EkB0hxy5N3LQpd5DLgkfVJDjqLIAxOd/yCHBjk0XZHoTNDnk0IIMhP68skB3yNu3i5sgvSlmBTtkjdZWcbICQuyCDInjEmS3/djYxx5GkHt7VQf97PDxGAG5w3hkwU0JPUCIr/0Z2ccQiIULsfWSzxUfTjrtSOP+9YQkTdw3Nqrx5rw1q6c1yFRkTc2nbMKGLKhBKMUht7bK9921S5wrXoJcFj2RA9RBm8rcIpVPJkDFC/LIiHryUyExYY+0bIKELdyENw6HTGkZLRYJE0Om7hgAtqINb0GzVG6z9AzxYxv9phAUvbGQaZs6ZMtKyCEnBHXIR8r7TawxZBoysoXEssTNwCumrAuyfq7/9KcxOOQsCVgHrTtkID9hi4oXZP0ktIfdXoSttKAXF3VvcQuym0N2iyFblllrn8F89Qm3WXoBBdnvuPXjsqGflVvIgopTXZ34/jwFOePl3qkgH3us3KbnwvCw6jabmsIJ8iw5ebLokG+7TYQuZs0y36T097Rv6p/4hHzu2mvLQ5Dffhv4yleAF17QfmGog966Zyx2/YtMKg4MmK9pFuSUCBOusAlbRUAFmf7fOJJ6FNshUzHs61Pfv75e/Ny3T5ybvoLskQWnwqgLcqSQBdwdspsgm/pYuApyxsu979kjRyS1te4OWTcJNTXhBHnmTLltC/IPfyh+btsm1h00HZtNU5N4T0CsTGPvb/t2WQM+frz4vOloLC+CfMklwNe/Dnz4w1q8XTvZ/4T3YgbWYuo7L2DZaAubNWvMVTAsyCkRRZCpOISNIdPhZBoOeeNGWeHR3OyMxRoFecypgbLgSYQs3BwyFWf6mlCCnPFMMToU7upSvwt6Luh9LIBwST1dkHt7gVdeMb+X6f/Tm+nYscB3vuN8vT1FP28OeedO4PnnxfbmzVpIThvp3YVLsQ+N2IOm4g3LbZo0C3JC7NsH3H+/dA5RBNlNHNxw68YWdtUQtzihjSmGTKcOt7Y6j90kyH8dmol3N/r3BXBL6gEeIQufkEGQGHJkQc54phgNV8yerZ5v9BzRO70BpTnkJUvU86zvH/+f43N3E2QA+Ju/Ad7/fvW5vAoy7UEBaMk4rQ56NbqK248/Ln7Sm2YNUT9O6iXEJz4BnHsuMH++iBelIcjU8bS2SrEKu2oI3Q/t12BjClnQKaAtLc5jpw7CnnUFAM8+6388Xg7ZGLLwCRkMDEhhqK1Vvw9bmJT9wV+Qd+wgYpTxTLGggqwn9IDSBPmZW1WV6ts26AjVeAlyoSBWOaf5lbwKst4DWxFSrQ56Td1BxV+99pr4rKggH3GE3GaHnBCPPip+rlgh7qZRBNlrGq8JfQjaIltnhIoj04uGLiZK9w2ogkxj3CZBpubwgx+U28884388oUMWPiEDva6Zlm9RkaCfg0mQ6+rk+ytVNBnPFKMlb3EK8siIepM68EC5vWkT8PSDak1dH8Y7QjVeggwAhx4qVju3sUU/74LsKPsj/Vh6xqkX0eLFashiPkmlsCD7MDQEPPSQc965F5YFbN8mx25vLPg/6H34ueLjNBzy+PHuF6If9KLp7FTjikD4kMXmzWqviYsvlttPP+1/PFtWyHhA+2UfVhyXUZB9QgZu8WN9f36CDLiELTKeKUYd8sEHxyfIenVEe7scbu/YATy372jl9X32ggzk+9B7IZv45jdFwuzss4HPjZaL560O2dMhE3budJqhxx9XHfLJJ8ttFmQfrr8eWLAAOOyw4I3id/33/2B4RP5Jb2ydit477i4+LidBHj9erdgAzA6ZlrzpDvmNN+RqCG1tapxw6VJvxzNw5y+wqkeOYadseFkZBhsF1CdkQKsn9Jh0EIdMb0SuceSUl9yhvPWW3C7FIetJPVNVBh2x7IU6KigKMvk+/ByyfSzd3aJxu/1f8+SQBweB119Xn3ObGGMqbXv0UdXg6Q45D2s05laQ//d/xc+9ew1LvrskjrZ99XvKy97Ae9A7KG1mUiELPUkTlyDTJaUAcwxZ/z0VZNpIqKNDXMT2NNvBQbXHhc5jX34IfaNd8rqwCtOxXhkGG2PIPiEDenPTBVkvr7KbLIVyyBmyfbsUh4YGYc6bmmRYZs8eWVLmV2WhO2STmE6Z4n4sfRjvCNXovZCDkidBXr7c2XHOzSGbBHnlSvnZtrUJw0PXFQxTT58UuRTkvXtRrBsEtLCFR+Jo+zp1TPU6DlNbbwYUZP1i91vBOa4Ysj6jMIhDpughC5qRth0PdQVeceTfvCvblF6A38jJ1qPDYGPIwidkQB0yPU5A3FvpPu3hca4FmRiDtw79aPHpgw8WCbJCQT3n7O+X3uRbW8VPr5BFaEEeO8URqgnikE3kqQ7ZtKhtGEGmzJolvp9p0+RzeQhb5FKQX39dFUFlqqNH4mj7VLU0oQedWAdZPxZUkMeMkRdKkBWckwhZNDe7O2TdhNLfU6GjF5BJkJcsMe9naAi4p0YKzAX4jWNHrmVvHiEDL4cMmMMWeh8Lm8wFWTMGyzfKuzCtkDGdC/TCt6sZYhXko9/nCNVEFeQ8OWSTIAcJWeg3fwA4aLQAo+oF+fnnxRRN03xyG71nr+KQPRJH2y++0vH0ksKJxe2gggyEC1ukLchBQxYUW5APO0w+5/ZRPvkksHVEfADTsQ7H40/iF2QYHGWmnpdD1vdpfxZBHHKknsilohkD2h8kjCDbghCrIBuaVFWqIAdxyJde6vz9rD2vAl1dmPbcb4vPVZ0g79snpjt+4xvAZZe5v85TkD0SR9sPf5/j6dVWV3E7jCCHSewlEUOOI2RBsT82Uw9cnd8QQ3z++EdQU4Aj/BBlpp5XUg+ILsjUIW/ZIlYSDrPai45liXPOc2Sk3c2WQ6owbYdpOhf0fhGAU5BpgskkpnSxW/09K1GQLcssyL295pmN9Ps/7zw54rWZ9fvvA2vWYBqkCm+8X2+OkT6pCvLatfKifP5596ymLshr1xLX4JE48hu60tiuH2EEWY8h04swah1yUg6ZDtHsXsaUkRHgblmYggvuvcwYfojSXChKyMJNkN1acF5wAXD55aKixK4wCcLevSKRfPnlYpbkrFkiFvzGGy7/QTMGVJCjOOSaGnXiDk1eBXHICxbI7UoU5LVr5UhowgTVWJjCFlSQZ84ETjtN/f1BA+KLnQp5d9x4z59iOtropCrI9MIZGDAPNYeHndMjLYusU+aROPIT5LRCFlT44yx7s0ukvGLIdMFQii3udXXyYrYs55DvuedkmWF7O3CKyxKEpTpk043DJPJhHPLAAPDUU2LbnhgUhKEh4MQTgY99TLhr++/fvFnU5ZpmzlFjMIgxWAE5K4y6VToD0eSQ6Q3SLWwRRJDppJ+kBDnLOmTaq+Ooo9QRgn4O79snn6utFTfYM89UXzMLImaqOORelwsrRTITZMA8ZP7LX8x3YiVs4ZI4ilOQgzrk/n7pxOrqxEUVJWShr8vX1CRuCiecIB6ffrqc2trYaBbd1lbxGn14VlenXvheYQsarvjIR9zblfrFkE2uOU6HbBLkdevUUVeQ6eGAqDahFzzl1VeBf/1Xwy+IMViFmRiCsLczZqh/h34u7N4t/7b6evXmHVWQx4wRN0571ZbBwWBx6CDkxSHTcIUuyLpDpu0C9t9ffC5nnCGfG1/owxSI/0QFeVMj6bubEakKsu40TYLstghnkBl7WQiyqVFMFEGmC5zW14t/hYIo0r/vPvHPplAwu2T74tbd54wZaiMV0+KWgBAzKsgXXOB+vF4O+e//XnwG116rPp9UUs/+3vVm/EEF+YEH5PZZZ4lRwk03yee+8x05JV9h1Bgs/52cEaIvp6SfC7o7pjfWMII8Z4483z7wAXE+UDeuu+RKE2R6Q9IdMg1X2L3N584FLrpIXAdfPH8lCqMXkOKQpx0V92GHJlOHbJqBRwWZXvReVRlu+9ehJ6wfQUMWpiL/KILsdsG0tADnnOOMG4cRZD0P6uaQly6VuarWVuHK3XCLIW/fLlaesCzgu9+VjpWukFwoqIJqY3LdYQRZrxoJKsi//73cvuYaMSr53OeAD31IPv/xj7vfmJcvl9t6Uyj9XDCVvNmEEeSWFnGT+Pa3gdtvF88lIch5qUP2cshBBLlQAH75S/EdXPvrI4ujm2k0hjzgkoBJkdyFLKggf/jDcjusQ+7qUn/X2KgmTfwI6pBN02CjTAwJe8HoAl1TI/+fnyC7OeSlS+X2ggWy2b0Jt5AFvXnZ67Tpz9MFU932aXLIflOndUHu6VF7eZhYs0Ym7hoa5E2oUAB+/GM5RXnDBuDf/s28j6CC3NdnTujZuE2fdjs33vte4MtfljfYSnXIvb3SkI0ZIxoheYUsTIJsU7xuRkc3U/euLv5u06Zw7XKTIFeCbFmqINMhc1hBpo1DgHDhCkB1yGEFOU6H7IYuyC0tcvgb1SGvWiW3Td3m3N5/927phPUbkL1/v3AFED1ksWPbMEYKtej5d+fKIM8953hKgYYrTj9d/bumTgV+8AP5+MknzfugPSz8HLKp5M3Gbfp00HOjUgWZJmcPPVTcuMKGLNwYO1bmXIaGhHFYtkxU6Xz+8+n3t8iVIL/zjnRSra1qIH7lSv8PJ05BpqLhFbKIK4ZsEnYvTIJsE9UhFytZ4Bxh6NTVyRHH8LAUEDdB9kvoAeEEua4OaGoQ2dQR1KIPzejZ5YyD+IUtaLiiGKIgU6JP+Pzxxd+bzgPLUqf5+8WQgzpkFmQJTbjay2IFdchBWmHT72HDBmGeH3sM+N73RMvONMlVUo+646OOEsJin9C7d3vNWiOLAAAeEklEQVQveW5ZyQlyUIdsv4d+YQQZBoW9YPQYMq2siMMh0567bpgSe/pkiqQcMgBMHJINm7djInrgvPq8BHnvXrmSBDAaItOmRE9aJ+3Ztm1OU7BlizzvmpqcK734JfUoSQmyZcXXXCiLjmi0Fvzww8XPsDFkL+j3cPvtqg7RMF4a5CqppwtyoaCusOsVtti1S1YpjBsnho50CJhmyKK2Vl44+sXgRhwhC5s0HDJgjiMHCVm4OWRTotCtlwUATByWO92OiXgH2kwaAH/+s7uz++Mf5e8OOWS0Kbs2JbpxdFU2QNQ567W4NH58yCFqNQvgk9T7hwuVboVugqw3nXLDTZBNFTxBGTNGltONjDi7raUBbblpT/13C1kMD6t5g7AO+b/+S/3da68FP844yFyQqXukgnz0aM9tulyNlyDTfU+cKESRxkHDCnJLi0w87drlfiK6hRrChi2SFGR9xp/uLoaG1BFIXZ0q2kGOwT7+tEIWADCxQf5yNbrQC/EhNBb2FmO5Q0PubUZp/LiYQDY0+JgEObTTb85e8WPAIMivy5vINGxQuhWGTerpuAlykOb0XmQZtrAsVZDf8x7xk/aD3rJF3nDWr5ezT9vb3SdRUej1oHd21PsvJ02mgjw4qJ7gfoLsVfqmCzIgvzwgvCAXCsFK3/IiyG4hi4kTnTHp+np5QluWEGI97uY2IcTtGMI45NhCFofJu8YrOLK43TFtQOlqZwpbWJZL/Nhgqagg6+eB7pB1HCGL1fIPKk7bHe1WGCRk4ZVfcBPkqOEKmywFeeNGeW3THuF1dfL6tCx5foUNVwDO0BHlzTed7QWSJFNBBlQHZRf2NzRIt6GELL5yh3ElY33ftiDb8SbAXQS8MIUtBgaAX/9axpbc1uxLWpD1Oz91yFRTZs82/389jkzjx0HCFUC4kEVYh+xXhwwAEw+X1l8R5MNacNJJ8nUmQX7rLfk3NzeTKeKGXimTamRg3EuQ/Rzyzp3AxmH5x9NJCejpSazKImpCzybLWmQ9XEEn0pgSe3SAE1WQGxqkhvT3B5sDERepCfLgoHmOvX3B6h+8ncGfueqx4vMrMdOxkrENvVDsD/Oyy4TDnjpVzB4Li6nS4tvfBi68UPQ+ePNNd4ccthbZGCd0WRkF8A5ZzJ4NfPGLYlXdb37T/H40JLF+vRo/DpLQ04/BTZDtGHUiDpkUVbzcJfs3d3TAIch6MkqfnVeMqxp6pUw6pqv4Wl2Q/UIWei+LfRB/RCP2oBlEKTs6jA45TEIuKUHO0iGbwhU2psReHA75058WNd42acaRUxNkt1aGtiDTTCr94Gf+7GvF7ZUYjV9oK+oCZoc8ZYpoMvPOOyJJGBZTpYV9IQ8OAj//eTIhi/Hj4bkyCuAdsgCAG28U5UJnnWV+P+qQN2yI5pCDxpAtK1oM2bKCCzK9oRxwgOjUZo9wtmxRVxsG1FIqutYgAEevlMlHySEHFeT+fpnXKBTMo5ExY8xxzGnYKFdhGe1WaBLkfftknqWhwXtyUyUKMtUF2ssbMCf2wpa8Aao5qa8H/vmf1fdKM46cmiC7TWu2HdSbb8rnDj1UbnesX4JaiCDOOszAXoyO6+yxyaiL3H7Fl4r/h16ohUK4GXoUPWSh15zef3+CMWSPlVEA75BFEOJwyEFCFv394rsPUmWhC/LgoBSj2lrn96ivXG3T0SEGFSfKtQkcy1XRv/egg+CJW8XNihXy+Do7nTcMG1P+Ymr9Dke3QlNSL4yYVqIgmyosbEwhiygO+fDDgVNPFdvf/KYoXaThzjQd8pi03shNkG2HTAWZOuS6zunoWNODVaPueDW6MBfLxVVnu8g9e7Ad0iJOfOcVgMQUo6KHLDZtUgXnlVfMdcj6diRB9lgZBfAOWQRBL30r1SG7CTIgvuMoIQsvdwyY+2EA0hkdf7xM3Ond3MJcuG7JXVr14xarB8S5oK9GMe1DxwB3qwXqJodczYI8MuI+cgb8HXJQQS4URD36nj3y86lah2wSZOqQcd11mFkjP+WVmCmXEiIucjvk1Tnx8d8iDvSQBXXHxeMhF6WbQw4SQ3ZcNB4rowD+IQs/9KReqTFkt5AFIO4hfo2FADFctGteaTMiwDzs9xPkgw+Wz9EbztBQuFpVN0GmIutVJmhyyKbMvimpF1WQ6f+LU5BD90T2yIP40dMjj33yZOdUc90hDw1FE2RA7QUDCA2yE4grVqR3I8pEkGnj9fXrhdjZd7jGRu2DXLgQs06V/+HticfJpYSIi1QEeXuAxhcB0IeqJkGmuCX1gjhkR1LPY2UUIF6H/NZbUmgaGpwnvhtBHfIbb8ikWmurFF0T9KLYLCfihXLI9vnlVsO+bp2sN5061T3UYOMmyF59KSjGkIXh9RXnkH3yIH7o8WO9B7ie1Fu8WN4wpk93Pz+C0NgoQ1kjI/7XflykJsj0RKZDj/XrVXc8Z46zBnbm2XI8+NJ5X8XwxaNLCRFrowhyezyRGD1kEUaQSw5ZeKyMApQeQ6YOmSa8bDMThCAxZECNwbnFj037pGGOoII8ZYp8rS7I9k0hrIty62tSiiCbHHLFCbJbHuTSSwO5Za/4MeAMWfzyl/LxhReaF3EIQxZx5EwcMg1JbNyo/rFKuGIUemHdcYcIul99NbDsqu8XlUkR5CsvjuWY9ZAFrTk1fdluMeQgq9kai/9dVkYBSg9ZuBXDB40f68dgJ+FMFyzt1hVGkKM4ZBp+mDRJ7YVi7y/sFHG3pF6SDjlPSb3IdchueRAgkFv2KnkD1M9w7VrgtyRSeXEMEpBFHDkTQZ42TZ7kw8NqW0PTB3/yyerJumkTcMstwHFfPxfbvnMH0NmpCvLFZ8dyzF4hi/POU19bU6OeuHTa9oMPAn/9q/d7lTx1+oFfhIrV1debxTFo/Fg/ht273WPldATkN0GH/u20mVRQQabTxN16oYR1yEmELNJwyPaIIDOH7BecN5SvUrxK3gCnQ7ZzDp2dIqFbKtQhV7QgT5qkxjAfk3M/jA55v/1EP4Krr1ZP/N27gT91XAisXo3tbTKDU0rsiELFo6dHxB4BUX71+c+rr21uVl3zEUfI+taREeAb3/B+r1IEuX7MMMb+wxWhY3WmRFQYh6yHLNwEma7+HKdDrqtz3ph0DTDFkakgh3XItONb3mLIY8bIxODIiBTPwPtwScBFFmRTHkTHxUUPD7tXXtk0NZknylx8cenhCkC9CVR0yGLiRFUMaKzQJMiA+HBuvlmI4iWXyOeXL3e23oxLkOmFSOPAs2cL105/b+ox8DU5pwXd3eqsLoo+GyuIINP3m2Rt9axZdoPGkW3icshuPRfCOGQ/QQac37UiyN3dmPnIrcWHK+8WtW80ZBHEITc2SrHs75eiFLcgl1plAZjDFoH24ZGAiyzINA/ihouLfvtt+Rnst597zbnpc4wjXAGIpJ79va9b579EXBxkktTTBdmmoUF1NCZqa4F58+Tj5cudrTep0ygFt33NnSvcCF163XTRnXQScPZo9GRkBPj6183v098vG5jU1QVrj7j//mJxSwC4fPh284u8Yngo0SF3d6Pp784vPty1eosiyG4rjsTpkAEPQR4VmJm9cjG2lXe/AnR3h3bIpkZT/f1yiFxT432j0c+N5mazsyvVIQMlCLLHRKSS6pDtPMhdd3lWDen4hStsaNgCEA2ejix9CgIAcY3T8ziNsEU2DvnMYzD9tzc7XmOqsDBBewYsW5aMOwbEhWi60Oz3P+cc+ZxbUo265F/8Qh2G2USN8f3hD6JK5brO28wv8InhRXbIo2LXtEl2Xdm9chN2/O6J4uNp08zfhZ8gU6EqSZBHBWYmZL3byuEOjPzrv0VqQKMneGl8u73d+7zVBdktoZqUINMEn2svDI+JSCZBvv56YQr0/sGu+FQN6fhVWNjoDjmucIVN2nHk9AR5nWzKOhHbsF/vcsdr3MIVOvSutXx5coIMmIdK9vufd574wmpqRCMjE8cfL1s7WpYq0DZB2yvqFAqjoupTs+yG7pDHjfMXTABFsaPNcXZbjdh5133Fxy0tZgeeWshiVGAUQcZMbOwZKPa2njw5+A1Qd8hBwxVANEGOUmUBmAWZrgrj2mrSYyKSPjGkrw+49lqx3y9/OUR7yoULccdXV+Oyj4/grYdWu4ox4F9hYaN/9hddFPBYApJ2HDmbpB62YTqcS04HFeSODhlv27xZbY8XtyCbBMQW5MZGsTz5pk0i3OYGFeFf/Uq2GbUpNQse1n3Y6A65qyugu7Cnb0PeZHehGTu3yivTTZDjrEMG1O+7vp4MYUcFpgM9qIGIZ63FDLw19X3F14eZyRWnILu9PgmHPDysCrK+xFQRj5u67pBXrJAhwj17nI2b3HjzTeDyy4E77xQJei9oRZOXINOQxZFH+i/OG5aKdMgDA8BuS4yVajGE8egrSZBratRm4HRl4TQEmb53TY2/yMybJ5eXB4BHHlF/X7IgA541y27oghk4oWdP3yaCvBtN2NkiXVYcDpnGpN2S9fT7PuAAMqllVGDqMYgDIO6AFmrw5Px/Kb4+TEVJGg45iaTehg1SPNvb1fdQ8LipmwSZElSofv97WaHy1FNqz2cdurSY143zhBPk9hVXBDuOMFCHbGofHDepCDJ1x63YgQJgFGSvO6EOvRPSBuRJhyw6OqKtukATgLogB10zLW5MDjkQo2I3DjIJtAdN2HGi/CNbWswx6jAOmRLEIStLVRGBoWGLx7ccUdyOyyF7rTgBZOuQ6WhMX8rLgctNXZ8Yogty0KE8Pe8HBsR6hyboSkKFgrpck86HPiRGnXfeCVx1VbDjCMOMGUJfdu5UVzRKitQFeSLEA2W1BIjqAlrE7wdN7NE105J2yFGHRLQv8aOPqmsJxuKQI6ALSWCHPCp2tZ0HFBcABYANrfLDcXPIbuVLNmEFmQq8Q2BHBWbmJ88oPrVkifx1GIesJ/WSjiFnIsgu+DrkB3uU+uWvnv8q5s4VK+vY7NsnXDHFbUXwzZulk25r8+59UiiIadIf/3jwKf9hKBREG9ewS8BFJX1Brt8NFAqo75yOtvFyJcdDDvH+4HWoINMFSPMqyEccIe/0W7ao04mzEuSGBlXQwgiULXZNbdI+0VilSZD9GgsB4QX5/PNlg6BPftL8GlpKSc+VPMeQ40rqpSHIr73YX6xffntNLb529xFYvlzkVey/4+mn1YVbATXUSAnz2VYaqQvypNOPKg6Jph8oA1pB48c2bsKYdMgiqiDX1ABnnikf0+Fb1CqLOKCumLarDAoN3/gJcpAKjrCCPHWqyDG++y5ZF0/DrbY9LUFuaFBryyvNIa+wZhUXjngK8kvYvh24916xrYfpALFogL60FqD2fvELB1Ua6TtkIpj0gg0ryLNnmysCknbIpnXTgkLDFm6CnKZDBkT50oEHAp/9rJpRDgo9XntqORBdkN3i814tMuvrvT83N0FOK6kHqHXqWQkybXsbBvrZb92q3ngBkSx9E+ICfhonK7/78Y/FT5Mgb9yoTmO3YYecMG6CTNe5O1n9Hn1pbDTHPPMasgBUQX7qqejD0jg55xzR4+EHP4j2/6mA0llcLS1O4Qmy8ndYhxwEkyBPmBCuQ14gQfZoxm6HU84+292pJlFlEbdDpqWIlNchyhGewXzl+YcfFsk7OyFWW6suIGoKW4RJmFYaqQiyaUVoAPjSl8S/730POOMM5//zw+RWkxTkyZO9M75+HHCALJnbt0/E1YDsqiziwM3RtrQ4Y9RJhCyCMHmyM44bKl4O9TzYtEmtAmhrg28z9uuvF87ywQfda73pmoFDQ+IcsWPetbXBWgLQz6+vT10ZJQ5BduM1HI7NaMNyqI7FstTY/vHHAwsWyMemxB4NWbBDTgA3hzxpEnDDDcDnPhdtumMagjxnjpwb7zYbLwx6tQUQ0SGXsDROnLgdr90wn4YtsnLIhYLTJYeJHwOqQ6bTpotVAD6L0gKiDNDrPC8UVNGltbh6N0E3qEPetk2KW6HgMSnEB7fPnn63r489Ds8Sd0xj5i/LdiI46yzR48WGHbJK+kk9n7KnMJjCB3ELck0N8PzzYpbRjTeWvj+65LwdVwstyCUujePYVwnCbnLItbXyeXrRZuWQAacgh3XI48aZmz4VHZzPorRBoYJMV8qmk5G8oIL81lsyaTZtWvTV1+vqzNUxtCf4a5NOxTNfuqf4+DOfMYeEzjpL7VX88stytRkbdsgJk1SviTQcMiAuxLlz42lactppshHNSy+JmNyu5XJcOf6fPuUvigHcWCBiEHaTIE+YID8rGucPklRym5EXtyCHdch6xzebomD4LEobFCrIixfL7fnzna81QQWZrm4SNVxhY/r8FyyQx7t+PXCfbGWCM890ThYdP17Ej1tb5SSw4WHghRfU13FSL2HSEuQ4W28mRUuLdAiWBTz8lWew6xXZjKN580p/UYzJjcUh7CZBps7ommtE9cYZZwAXXOC/v9pasyj79Tn3o1SHDPgIcsQGTzo0sUcFmQ7zvXArm4xaYWFjEuRDDlGro+gSZyedJPpWUE4/Xbp0R9iCjNQ2LpeCwSGLBHBL6pVKW5s6DE7CHScRq6Vx5Kt+cDjeGJF3lmbs8hfFmNxYHMJuCjHQBVfnzBGTYB57LLiomvaZtUMGfAQ5YoMnHWooVq2S26UKctwO2b4cTK0x58wRye+jj1Z7E9Pz/sQT5faz/7u2OFIbsMZg24i4kGsKI8G6D1YQZe2QAdUlxy7IccZqCZ/4hBSdndYEvAs5Liu2tPQSxZjcWBzCbnLIYVfA1klDkKM4ZFNSUhlSR2jwpGMa4XV2mqehm3CLwcctyB0d4lhNtet2eKVQAL71LeH6Dz5YLDZtozjkV8bBGh2pvQvZvq29sDVQf/RKoqyTekDCghxXrFajs1Ms7GoaRhYF2UsUY3JjcQh7EoJs2mepgtzZKY9r2rRgFR86ng45JkyCHNQdA8K5mj6/uAX5oIPET5NDpnMKFiwQ1/+yZWoo6+CD5ee5dWQS/gIxTXQTMSdTR5wNyCqdxAWZrkE2Zky0Tmle0EqL2AU5rlitgaOPFtUbx81UK+1bsDOYKMbgxuIQdr+QRRT0fdbURK8QsKmvB+64Azj3XOAnP4mWoC0HQQbMYYukBNnLIduMHets/FMoqBNEXoPYERXkaWNdVs2tYBIXZD1cEefyKoCIS9n7DHvi+hJXrNaF/fYD/vhaG/7+lLdRhwF8Cj9Ca2drNLcblRKFPY2QRWNjPOfNRz4ieivQiQlhYEGWzJ4tfu6/v/p9T5kixdoPJSFYJwR5I2QWb+q8EjORZUjigpxUQs/m8MOBJ54AfvYz4AtfiHnnccVqPRg3DrjjyVno21eP26xPRXe7GZGWIOeBNAR57DZ1mN7UMIgjjnB5sQu6INfUlF6toF8GB33rCqCmBoUDu3D4VDlT5uSTg988aW318uM+DnR2YhMV5BN8VjyuQFJ3yAoxVTCccopIGMRe8hZXrDYAeS/Xc6OaBNkUd9ZXPS6J7m40vLFUeer4wWcw5n/CXRe6IE+fHq61rQlHyGLzs8VE99ErZePj970PgaH5n+X9BwKrV2PjNdcXn6u2kjcgZUFWHEZCFQyxE0estoJJI4acF0HWHfKkSaXHthUWLULDsJpEPmnkqdBJZF2QSw1XAOp3UMCIsgrLPw3dgBMbluKDHwy3jBIVZHtWYSaTQnLShgDI0iEnVMHApEtsDplcFM3dtyq/yqsgxy4YPT1ogLrQ3El4NnQSOWlBnoG1GEuOswtr8OzAcXjggXBJ+/Z2WXmxa5eY7Zd6L+ScGcPsBDnBCgYmPWIRZO2iaO5V46hVI8gdHQ5BPgFLQieRkxbkg2BYZjpCortQ0MIWyzNwyDkzhtkJcsIVDEw6xBKy0C4Kupo1UPq06bhIXJCvuw5jxwwVHx6KNzBx3EDoJHLigly7Wv1lCYluPWyRukPOmTHMrsoihQoGJnliccjayV+cHDNKXhxyc7MaM45dkBcuRMMHTi0+PKn51UhJZF2QS+1jAahToN/3qUNiS3RTQX7lFWDHDrFdWxtt8k5ocmYMExfkCRNEaLClRfuAU6xgCE2Ogvx5Z+xYZ5lTaEHWTv68CrLe8S2JIfX8K2Rx7kV3/22k6yEJh3z++cCPfgT88IfA3948P7ZENy19e/JJud3enswq0g7yZgwtywr879hjj7Uqnrvusqxx4yxLhPjFv3HjxPOMkeZm9ePatSvkDrTP/B78jbK/v/u7RA47EnPnyuO6/fb49z8yYlmLF1vWM89E38dtt6nfx/r1sR1e7Cxbph6r/e+oo1I8iLvusqzOTssqFMTPBK51AC9aATQ2lV4WZUXOgvzlAI0ju7XP9EQbLTVPUeMgeXHIQPIOuVAQPbNLmXU6/pWni9t1GMDUR/M7wps1y1wjnWof5ByVtrIg6+QsyF8O0DhyS0vEac7komi+9+fKr/IkyHZP5/Z24NRTvV+bCd3dGH/7d4oP98c61Hw2h/X9o9TVCVHWqcZJIQALspOcBfnLAV2QSyWvE0MA4B//EXjtNWDFCvfew5myaBGm98tGyrPx19yP8EzLUwV2yBWW72FB1slbkL8MaN67ubjdsvaNki+KPAsyIFpO6qtYJ0IUsenpwVF4Gf+A7+MYLMVX8dXi83nFtBRbIIecs0kdccCCrJPn6o880t2NprdfLT5sGdxc8kWhl9LlTZAjEVZco4pNRwcKAL6Pa7AU83ASnis+n1dMghzIIVdgvocF2UQSQf4KG1oVWbQITSO9xYct2FnyRZF3hxyaKOIaVWzKcIRnCll4OmT7Wlqzxvz7HI8G/GBBToMKHFoV6elRZta1YGfx+ag0NEBZuqfsBTmKuEZNLpfhCM8YQ37hfrOBodeSGzkeDfjBgpwGFTi0KtLRoUzkKApyCRdFoQA0NwwUHzf+32vK++YVRVxLSS7nqIwrCJMni6oVytSvXWU2MKZriZLz0YAfLMhpUMmldNddh6YxUjwDL0HlRXe3kigct6WnvEcUUcS1DEMPpUDjyGMwiEl716ovsA2M1zVTBqMBP1iQk4LGjN3mgJbx0KrIwoU49tPHFh8e176m9Iti0SI0W33Fh43YW94jiijiWoahh1KgYYspeBc1sJwv6ulxv2Y6O8tiNOBHiesIMEbsOJc9tBoedr6mgtzOx26aj8H3iiL/cy/+Uem3+Z4eJQzSiL3F58sSWyRsh9fRIb57P/FYuLDsBSYo1CFPq98GDBheZH9u9NoCKupa4l4WSdDZaZ6gX1ub6Hz5iqGz0/oibrAAy5qAHVYvRptldHZmfWRMQjzyiLxMLjhujXc/mRR6T8QNAvayKIjXBmPevHnWiy++mNzdoVKoqRGnkU6hIBItjDfd3Rj41NV4aO8pOBKvoAPvCBdUwUP2aseygM98Bnj1VeCWW4BjlnWHH1HkmEKhsNSyrHl+r+MYchLw9OvSWLgQ9T/6L5zb+Ro6CmsrPn5aFiRRR0/2WTiwC7ed2o0lS4BjjkHZVYrERhAbbf/jkIWG29CJW3hWLmU4XC6ZJM7nKrtGEDBkwYIcFb8Tqhov3EqnykSkiFtOpJSYfhL7zDFBBZljyFFxm7ppl98w+aI7hphktX7nSeREqizPwjHkpKnkyR5xkoceHnFNXa/W7zyJnAjnWYywIEeFTyh/8tLDI66p69X6nScxa7DKZiIGJkhcw/7HMWRCtcYTw5CXOGGhYD6OQiHcfqr5O08iJ1JFeRZwUi8FquiEikRcQlgqcd4YyvE7L8djrjCCCjIn9ZjkyEsSTJ/KDlTPRJNq/ttzBCf1mOzJS5ywyhr1KFRy69cKJHlBzkOWncmGPAlhtc78qtbKkDIl2W5v+nDJzrID1XNBVDtV1LEsl3R0mMNGlV4ZUqYk65B5uMQw2ZKXsBETiGQFmYdLDJMteQobMb4kK8jVWkjvB8fVmTSp1vh5GZKsIPNwyUleZq8xDJM7khVkHi454bg6wzAuJF/2VmnDpVLDDRxXZxjGBZ4YEoY4wg0cV2cYxgUW5DDEEW7guDrDMC6wIIchjnADx9UZhnEh2Zl6lUZcs5549hrDMAbYIYeBww0MwyQIC3IYONzAMEyCcMgiLBxuYBgmIdghMwzD5AQWZIZhmJzAgswwDJMTWJAZhmFyAgsywzBMTgi16nShUNgMwDAzgmEYhvGg07Ksdr8XhRJkhmEYJjk4ZMEwDJMTWJAZhmFyAgsywzBMTmBBZhiGyQksyAzDMDmBBZlhGCYnsCAzDMPkBBZkhmGYnMCCzDAMkxP+PxT5xgo3fXuvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn import linear_model\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "\n",
    "\n",
    "regr = linear_model.LinearRegression()\n",
    "regr.fit(X_train, y_train)\n",
    "\n",
    "y_pred = regr.predict(X_test)\n",
    "\n",
    "# The coefficients\n",
    "print('Coefficients: \\n', regr.coef_)\n",
    "# The mean squared error\n",
    "print(\"Mean squared error: %.2f\"\n",
    "      % mean_squared_error(y_test, y_pred))\n",
    "# Explained variance score: 1 is perfect prediction\n",
    "print('Variance score: %.2f' % r2_score(y_test, y_pred))\n",
    "print(r2_score(y_test, y_pred))\n",
    "\n",
    "# Plot outputs\n",
    "plt.scatter(range(X_test.shape[0]), y_test, color='red')\n",
    "plt.plot(range(X_test.shape[0]), y_pred, color='blue', linewidth=3)\n",
    "\n",
    "plt.xticks(())\n",
    "plt.yticks(())\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cross Validating - Run 1 out of 5\n",
      "Fold 1 out of 5\n",
      "Fold 2 out of 5\n",
      "Fold 3 out of 5\n",
      "Fold 4 out of 5\n",
      "Fold 5 out of 5\n",
      "Cross Validating - Run 2 out of 5\n",
      "Fold 1 out of 5\n",
      "Fold 2 out of 5\n",
      "Fold 3 out of 5\n",
      "Fold 4 out of 5\n",
      "Fold 5 out of 5\n",
      "Cross Validating - Run 3 out of 5\n",
      "Fold 1 out of 5\n",
      "Fold 2 out of 5\n",
      "Fold 3 out of 5\n",
      "Fold 4 out of 5\n",
      "Fold 5 out of 5\n",
      "Cross Validating - Run 4 out of 5\n",
      "Fold 1 out of 5\n",
      "Fold 2 out of 5\n",
      "Fold 3 out of 5\n",
      "Fold 4 out of 5\n",
      "Fold 5 out of 5\n",
      "Cross Validating - Run 5 out of 5\n",
      "Fold 1 out of 5\n",
      "Fold 2 out of 5\n",
      "Fold 3 out of 5\n",
      "Fold 4 out of 5\n",
      "Fold 5 out of 5\n",
      "\n",
      "Overall R2: [-670.4095204  -671.45443974 -673.60734422 -671.40192667 -675.31208185]\n",
      "Mean: -672.437062575318\n",
      "Deviation: 1.77669784591671\n"
     ]
    }
   ],
   "source": [
    "import numpy as np \n",
    "import pandas as pd \n",
    "from sklearn.utils import shuffle\n",
    "from sklearn.model_selection import KFold\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "### 交叉验证\n",
    "def cross_validate(model, x, y, folds=5, repeats=5):\n",
    "    \n",
    "    ypred = np.zeros((len(y),repeats))\n",
    "    score = np.zeros(repeats)\n",
    "    for r in range(repeats):\n",
    "        i=0\n",
    "        print('Cross Validating - Run', str(r + 1), 'out of', str(repeats))\n",
    "        x,y = shuffle(x, y, random_state=r) #shuffle data before each repeat\n",
    "        kf = KFold(n_splits=folds,random_state=i+1000) #random split, different each time\n",
    "        for train_ind, test_ind in kf.split(x):\n",
    "            print('Fold', i+1, 'out of', folds)\n",
    "            xtrain,ytrain = x[train_ind,:],y[train_ind]\n",
    "            xtest,ytest = x[test_ind,:],y[test_ind]\n",
    "            model.fit(xtrain, ytrain)\n",
    "            #print(xtrain.shape, ytrain.shape, xtest.shape, ytest.shape)\n",
    "            ypred[test_ind]=model.predict(xtest)\n",
    "            i+=1\n",
    "        score[r] = R2(ypred[:,r],y)\n",
    "    print('\\nOverall R2:',str(score))\n",
    "    print('Mean:',str(np.mean(score)))\n",
    "    print('Deviation:',str(np.std(score)))\n",
    "    pass\n",
    "\n",
    "cross_validate(regr, X, y, folds=5, repeats=5)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
